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INTRODUCTION 


Chemical  process  design  requires  accurate  values  of 
thermodynamic  properties  for  the  various  chemical  compounds 
involved.  Reactor  and  separator  design  require  activities 
and  fugacities;  calculation  of  process  heating  requirements 
calls  for  enthapies;  and  equipment  sizing  demands  knowledge 
of  the  volumetric  properties.  Due  to  increasingly  stiff 
competition  and  scarce  raw  materials,  the  ultimate  costs  of 
over-design  have  become  prohibitive.  Consequently,  the  need 
for  accurate  and  reliable  values  of  thermodynamic  properties 
has  intensified  in  recent  years. 

Ideally,  the  design  engineer  would  have  experimental 
data  available  for  the  substances  of  interest  at  all  the 
temperatures  and  pressures  of  process  operations.  However, 
due  to  the  vast  number  of  state  conditions  typically 
encountered  and  the  expense  of  establishing  and  operating  a 
physical  properties  laboratory,  the  amount  of  raw  experi- 
mental data  is  limited.  There  is  also  a need  for  continuous 
data  to  use  in  control  and  optimization  schemes.  Thus, 
despite  the  fact  that  reliable  experimental  data  are 
preferable  to  either  theoretical  predictions  or  empirical 


correlations,  the  design  engineer  is  often  forced  to 
estimate  key  thermodynamic  quantities. 

Although  several  promising  fundamental  theories  exist 
for  estimating  thermodynamic  properties,  at  present  the 

serves  as  an  impetus  in  the  development  of  physical  property 
correlation  and  estimation  techniques.  The  primary  objec- 
tive of  these  efforts  is  development  of  completely  general 
procedures  for  accurate  and  rapid  calculation  of  fugacities, 

multiphase  fluid  systems.  However,  limited  success  of  the 


although  the  ultimate  goal  is 


The  present  work  focuses  on  the  formulation  of  a "group 
contribution"  method  for  predicting  equilibrium  thermo- 


that  individual  segments  (i.e.  groups)  of  molecules  possess 
thermodynamic  properties.  The  groups  may  refer  to  organic 


radicals,  but  may  be  more  specifically  defined  to  be 
individual  atoms  or  interaction  sites.  Molecular  properties 
are  then  given  by  appropriate  summations  over  the 
constituent  group  properties.  The  attractiveness  of  this 
approach  rests  upon  the  observation  that  the  vast  majority 
of  compounds  of  commercial  interest  are  comprised  of  a 
relatively  small  set  of  distinct  functional  groups  (e.g. 
methylenes , carboxyls , amines  and  hydroxyls ) . In 
particular,  since  molecules  in  a homologous  series  consist 
of  the  same  groups  in  different  proportions,  data  for 
several  members  of  the  series  can  be  used  to  deduce  values 
for  group  thermodynamic  properties.  These  in  turn  can  be 
used  to  predict  properties  of  other  members  of  the  series. 
This  concept  is  readily  extended  to  mixtures.  While  the 
present  formalism  for  group  contributions  is  written  for 
(and  is  applicable  to)  mixtures,  its  consequences  will  be 
explored  only  for  one-component  systems. 

Chapter  II  surveys  the  methods  available  for  estimating 
thermodynamic  properties  of  fluids  from  molecular  theories, 
which  range  from  ab  initio  calculations  to  empirical 
correlations.  In  Chapter  III  a statistical  mechanical 
analysis  of  the  most  successful  and  widely  used  group 
contribution  method  (the  solution-of -groups  method)  is 
presented  in  order  to  expose  all  the  inherent  molecular- 
level  assumptions.  In  light  of  this  analysis,  chapter  IV 
outlines  the  basis  of  an  alternative.  Here  the  concentric 


shell  intermolecular  potential  is  introduced  and  its 
properties  explored  within  a group  contribution  framework. 
In  Chapter  V group  contribution  equations  of  state  are 
developed.  The  first  is  of  the  Van  der  Waals  form,  the 
second  is  based  on  the  statistical  mechanical  pressure 
equation  along  with  a mean  field  approximation.  Then  in 
Chapter  VI,  the  results  of  model  calculations  based  on  the 
proposed  group  contribution  equation  of  state  and  molecular 
dynamics  simulations  using  the  concentric  shell  potential 
model  are  reported.  Finally,  in  Chapter  VII,  the  primary 
conclusions  of  the  study  are  summarized  and  recommendations 
for  future  work  are  offered.  There  are  several  appendices 
which  contain  details  of  the  derivations  presented  in  the 


CHAPTER  II 

CALCULATION  OF  THERMODYNAMIC 
PROPERTIES  FROM  MOLECULAR  THEORIES 

Molecular  theories  for  calculating  bulk  thermodynamic 
properties  invariably  use  statistical  mechanics  to  link 
molecular  quantities  to  bulk  fluid  behavior.  The  link  to 
equilibrium  properties  for  most  substances  of  interest 
involves  various  formulations  of  classical  statistical 

of  either  the  molecular  partition  function  or  the  molecular 


tion  functions  (O'Connell,  1971;  Reed  and  Gubbins,  1973). 
Each  approach  has  been  developed  and  used  to  calculate 

requires  detailed  knowledge  of  the  intermolecular  forces. 
In  what  follows,  available  methods  are  critically  reviewed 


and  prohibitive  for  the  class  of  problems  encountered  in 
chemical  process  design.  Fortunately,  classical  mechanics 
along  with  knowledge  of  the  forces  seems  adequate  for 
estimating  properties  in  most  systems  (Gray  and  Gubbins, 
1984).  Exceptions  include  low  temperature  processes  and 
those  involving  the  light  elements  (e.g.  hydrogen,  helium 
and  neon).  In  such  cases,  quantum  corrections  can  be 
imposed  on  the  classical  results  (McQuarrie,  1976).  An 
outstanding  question  is  the  impact  on  bulk  properties  of 
vibrational  transitions  and  relaxation  occurring  at  higher 
temperatures  and  how  quantum  aspects  can  be  superimposed  in 
such  cases  (Berens  et  al.,  1983). 

The  classical  equilibrium  fluid  is  a conservative 
system,  thus  the  forces  are  obtained  directly  through  the 
gradient  of  the  potential  energy  (Goldstein,  1980).  Quantum 
mechanical  techniques  may  yield  potential  energy  hypersur- 
faces for  thermodynamic  property  calculations  via  classical 
prescriptions.  Established  computational  methods  such  as 
molecular  orbital  theory  and  the  so-called  ab  Initio  methods 
may  be  used  to  generate  tables  of  numbers  representing  the 
hypersurface  at  various  points  (Clark,  1985).  The  name  of 
the  latter  method  is  somewhat  of  a misnomer  since  neglect  of 
relativistic  effects  and  use  of  the  Born-Oppenheimer 
approximation  are  common.  These  numbers  can  then  be 
interpolated  (and  extrapolated  when  necessary)  by  least- 
squares  fitting  to  a function  having  characteristics  of 


this 


surfaces  and 


effects  due  to  many-body  interactions  on  the 
resultant  properties  (Egelstaff,  1987)  must  first  be 
resolved. 

A common  (and  practically  necessary)  assumption  is 
that  of  pairwise  additivity.  Here  the  system 
configurational  energy  is  assumed  to  be  the  sum  of  all 
two-body  interactions.  Potential  parameters  deduced  from 
experimental  data  using  the  pairwise  additive  assumption 
implicitly  account  for  many-body  interactions  and  can  be 
used  to  generate  an  "effective"  two-body  potential.  For  a 
pair  of  polar  anisotropic  molecules,  the  multipole 
expansion  provides  a means  for  determining  the 
electrostatic  interactions  which  may  then  be  added  to  the 
nonpolar  pairwise  additive  contributions  (Reed  and 
Gubbins,  1973;  Gray  and  Gubbins,  1984).  The  expansion 
is  based  on  the  idea  that  any  function  dependent  upon  the 
molecule's  charge  distribution  may  be  expressed  outside 
the  molecule  as  a Maclaurin  series  in  ideal  classical 
raultipole  moments  (dipole,  quadrupole,  etc. ) . The 
usefulness  of  the  series  is  limited  by  the  requirement 
that  molecular  center  of  masses  be  separated  by  distances 
greater  than  the  distance  from  any  molecular  charge  center 
to  its  respective  center  of  mass.  This  limitation  becomes 
severe  for  dense  fluids  since  at  liquid  densities  molecular 
charge  clouds  overlap  to  a degree  and  structure  is  deter- 
mined by  the  resultant  repulsive  forces  (Andersen,  Chandler 


and  Weeks,  1976).  Due  to  overlap,  the  charge  distribution 
is  distorted  and  is  no  longer  simply  a summation  of  isolated 
multipoles.  Use  of  quantum  calculations  for  short-range 
interactions  and  multipole  series  for  the  long-range  inter- 
actions is  possible,  but  "matching"  the  two  along  with  the 
difficulties  mentioned  above  has  hampered  implementation. 

Other  formal  methods  have  been  proposed  for  calculating 
potential  energies  as  well  as  key  statistical  mechanical 
quantities.  These  include  expansions  of  the  potential  in 
terms  of  spherical  harmonics  (Gray,  1968 j Gray  and  Van 
Kranendonk,  1966),  ellipsoidal  harmonics  (Stiles,  1979),  and 
Chebyshev  polynomials  (Monson  and  Steele,  1983)  as  well  as 
spherical  harmonic  expansions  of  the  pair  correlation 
function  (Street  and  Tildesley,  1977).  While  harmonic 
expansions  are  elegant,  possess  formal  rigor,  and  can  serve 
as  a convenient  means  of  obtaining  pair  correlation  func- 
tions using  computer  simulation  (Gray  and  Gubbins,  1984), 
their  practical  utility  is  yet  to  be  demonstrated. 

Thermodynamic  Properties  from  Computer  Simulation 

In  the  last  two  decades,  computer  simulations  of  dense 
fluids  have  become  widespread.  Traditionally,  their  use  has 
been  to  provide  essentially  exact  experimental  "data"  on 
well-defined  models  which  then  are  used  to  test  theoretical 
calculations  using  the  same  potential.  Removal  of  uncer- 
tainty in  the  interparticle  potential  allows  theoretical 


l ambiguously  evaluated. 


still  (implicitly) 
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dependent  upon  molecular  orientation.  The  individual  sites 
are  typically  assumed  to  interact  via  some  model  potential 
energy  function  le.g.  Lennard- Jones , Kihara,  n-6,  exp-6, 
Coulomb,  or  some  combination  of  these),  often  ab  initio  or 
molecular  mechanics  calculations  (Clark,  1985)  are  utilized 
to  determine  parameters  for  intramolecular  interactions  such 
as  rotational  potential  energy  barriers.  Once  intra- 
molecular and  site-site  intermolecular  potentials  are 
determined,  the  simulation  is  performed.  The  alternative  to 
the  site-site  potential  is  to  apply  one  of  the  potential 
models  to  the  entire  molecule. 

A number  of  fluids  and  mixtures  have  been  simulated 
(see  Quirke,  1986,  for  a partial  listing)  and  calculated 
properties  compare  favorably  with  real  experimental  data, 
such  calculations  usually  begin  by  using  spectroscopic, 
electron  structure,  and  molecular  mechanics  calculations  to 
fix  a priori  bond  lengths,  partial  charges  on  sites,  and 
rotational  potential  energy  barriers.  Then  a number  of 
simulations  are  performed  in  order  to  deduce  the  optimal 
values  for  the  remaining  parameters  (e.g.  site  Lennard- Jones 
size  and  energy  parameters ) . Optimal  values  are  found  by 
comparison  to  experimental  data,  since  actual  data  on  the 
system  of  interest  are  required,  these  methods  are  not  truly 
predictive.  However,  they  may  be  useful  in  interpolating 
between  available  data  points  or  to  extend  the  database  to 
include  temperatures  and  pressures  under  which  laboratory 
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ar  radial  distribution  function.  A postulate  of 
classical  statistical  mechanics  is  that  ensemble  averaged 
values  correspond  to  thermodynamic  quantities. 

with  the  pairwise  additivity  approximation,  expressions 
for  the  thermodynamic  energy  and  pressure  involve  only  the 
pair  correlation  functions.  These  take  the  form  of  the 
familiar  energy  and  pressure  equations  (Reed  and  Gubbins, 
1973).  Thus  knowledge  of  the  pair  correlation  function  and 
intermolecular  potential  are  sufficient  for  obtaining 
equilibrium  thermodynamic  properties.  Also  note  that  it  has 
been  proven  (Gray  and  Gubbins,  1984)  that  for  each  distinct 
pair  potential  function  (differing  by  more  than  a constant) 
there  is  a unique  pair  correlation  function,  that  is,  there 

A number  of  integral  equations  for  calculating  pair 
correlation  functions,  based  on  knowledge  of  the  pair 
potential,  have  been  proposed.  The  Kirkwood  equation 
(McQuarrie,  1976)  and  the  Born-Green-Yvon  equation  (Hill, 
1956)  consist  of  a hierarchy  of  nonlinear  integral  equations 

(m+l)th  molecular  correlation  function.  The  hierarchy  can 

a single  equation  involving  the  pair  potential  and  the  pair 
correlation  function.  This  superposition  approximation 
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example,  Che  triplet  correlation  function  to  be  written  in 
terms  of  products  of  pair  correlation  functions.  The 
resulting  integral  equations  prove  difficult  to  solve. 
Numerical  solutions  have  been  obtained  for  model  potential 

for  simple  fluids  with  spherically  symmetric  potentials  but 
can  be  generalized  to  include  anisotropic  molecular  fluids. 
Due  to  computational  difficulties  and  lack  of  analytically 
simple  and  accurate  pair  potentials,  relatively  little 
effort  is  currently  expended  in  calculating  thermodynamic 
properties  from  this  class  of  integral  equations. 

A convoluted  integral  equation  was  introduced  for 
spherical  molecules  by  Ornstein  and  Zernike  in  1914  (see 
Frisch  and  Lebowitz,  1964).  This  equation  defines  a new 
quantity,  the  direct  correlation  function  (DCF),  in  terms  of 
the  pair  correlation  function  and  has  led  to  several 

an  additional  "closure"  relationship  which  relates  the  DCF 
to  the  pair  correlation  function.  With  such  a relation,  the 

integral  equation  for  the  pair  correlation  function.  Useful 
closure  relations  include  the  PY  (Percus  and  Yevick,  1958)  j 
the  hypernetted  chain,  or  HNC  (Rushbrooke  and  Scoins, 


1953) 
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to  calculate  molecular  pair  correlation  functions.  In 
presently  available  ISM  models  it  is  assumed  that  the 

The  reference  ISM  (known  as  RISM)  utilizes  hard  sphere 
potentials  for  the  site-site  interactions  and  has  been 
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es  (see  Friedman,  1985)  require 
knowledge  of  the  pair  correlation  function  and  the  inter- 
molecular  potential  to  be  utilised  in  a truly  predictive 


In  perturbation  theories  of  fluids,  properties  of  the 

are  known,  by  an  expansion  in  powers  of  some  perturbation 

difference  between  the  intermolecular  potential  energy 
functions  of  the  real  system  and  that  of  the  reference 
system.  The  perturbation  terms  (first  order,  second  order, 
etc . ) are  then  given  in  terms  of  multidimensional  integrals 
involving  the  perturbation  parameter  and  various  order 
correlation  functions  of  the  reference  system.  Details  of 
the  derivations  are  omitted  but  can  be  found  in  reviews  of 
the  subject  (Hansen  and  McDonald,  1976;  Boublik,  1977; 
Wertheim,  1979;  Gray  and  Gubbins,  1984).  Full  accounts  of 
perturbation  theories  for  simple  fluids  (i.e.  atoms  with 


srically  symmetric 


anistotropic  contributions.  This  is  accomplished  by  using 
the  Mayer  £ function  as  the  expansion  functional  (Perram  and 
White,  1974j  Smith,  1974).  Applications  to  dipolar  and 
quadrupolar  hard  spheres,  dipolar  and  quadrupolar  L-J 
molecules,  homonuclear  diatomic  L-J  molecules, 
with  embedded  guadrupoles,  homo-  and  heteronucl 


diatomi 


dumbbells,  linear  triatomics  of  fused  hard  spheres,  hard 
spherocylinders,  and  four  center  tetrahedral  L-J  molecules 
have  been  reported  (see  Gray  and  Gubbins,  1984,  for  a 
review) . For  molecules  with  significant  anisotropy  such  as 
spherocylinders  or  dumbbells  with  large  elongations, 
convergence  is  slow  and  poor  results  are  obtained  for  the 
thermodynamic  properties  especially  at  high  densities. 

Several  other  perturbation  theories  based  on 
nonspherical  reference  potentials  have  been  proposed.  Blip 
function  theory  (Sung  and  chandler,  1972)  and  Bellemans' 
theory  (Bellemans,  1968)  use  nonspherical  repulsive  poten- 
tials and  properties  are  expansions  about  those  of  hard 
nonspherical  molecules.  Application  to  molecular  fluids  of 
arbitrary  shape  and  size  requires  detailed  knowledge  of  the 
geometry  of  fused  hard  spheres  so  that  computer  simulation 
may  be  employed  to  generate  properties  of  the  reference 
fluid.  This  proves  to  be  a nontrivial  geometric  problem, 
although  some  progress  has  been  made  for  molecules  comprised 
of  3,  4,  6,  and  12  fused  hard  spheres  (Lustlg,  1985,  1986). 
The  blip  function  theory  shows  agreement  with  simulation 
data  but  deteriorates  with  increasing  anisotropy  (street  and 
Tildesley,  1976).  Bellemans1  theory  proves  accurate  for 
thermodynamic  properties  but  is  poor  for  structural 
properties  when  applied  to  the  hard  dumbbell  fluid  (Kohler 
et  al.,  1979).  A recent  study  (Lustig,  1987)  of  CC14,  CF4, 
neo-C5Hi2  and  SF6  using  a fused  hard  sphere  reference  system 


manages  to  obtain  an  accurate  value  for  the  free  energy. 
However,  the  first  order  perturbation  term  for  the  free 
energy  is  larger  in  magnitude  than  and  of  opposite  sign  to 
the  reference  term.  This  raises  questions  of  convergence 
and  seems  to  violate  the  spirit  of  perturbation  theory: 
properties  are  expanded  about  a known  state,  in  powers  of 
some  small  parameter.  Other  reference  fluids  such  as  convex 
hard  bodies  (Boublik  and  Nezbeda,  1986)  have  also  been 
employed  with  limited  success. 

It  is  apparent  now  that  realistic  nonspherical  refer- 
ence potentials  are  needed  for  perturbation  theories  of 
molecular  fluids  in  order  to  insure  (and  speed)  convergence. 
Site-site  potentials  are  convenient  and  may  be  used  directly 
in  site-explicit  perturbation  theories  (Naumann  and  Lippert, 
1981:  Naumann,  1982:  Nezbeda  et  al.,  1983).  Alternatively, 
they  may  be  used  in  integral  equation  theories  (e.g.  RISK) 
or  in  computer  simulations  to  generate  properties  of  the 
reference  fluid.  This  increases  computational  effort  and 
costs,  and  decreases  final  numerical  accuracy  (due  to 
truncation  error,  if  nothing  else).  Thus,  the  fundamental 
difficulties  facing  widespread  use  of  perturbation  theory 
for  calculation  of  thermodynamics  are  related  to  the  lack  of 
simple,  yet  accurate,  reference  potentials  and  computational 
costs  becoming  prohibitive  when  the  available  "realistic” 
reference  potentials  are  employed. 
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leads  to  so-called  n-fluid  van  der  Waals  theories  intended 
for  simple  mixtures.  These  have  been  tested  for  hard 
spheres  and  L-J  mixtures  (Henderson  and  Leonard  1970;  1971) 
and  are  accurate  provided  the  size  differences  of  the 
species  present  in  solution  are  not  large.  The  1-fluid 
model  is  the  most  accurate.  Mixtures  with  components 
differing  greatly  in  size  require  additional  statements, 
such  as  the  "mean  density"  approximation  (Mansoori  and 
Leland,  1972).  However,  these  prove  to  be  of  limited 
accuracy  for  large  size  ratios  and  dilute  solutions  (Ely, 
1986a) . Success  of  conformal  solution  theory  for 
anisotropic  molecules  is  unlikely  since  the  assumption  that 
the  intermolecular  forces  are  conformal  for  all  species  is 
unlikely  to  hold  since  it  assumes  that  the  dependence  on 
reduced  (or  scaled)  distance,  energy,  and  orientational 
parameters  is  described  by  a common  functional  form. 

Macroscopic  corresponding  states  correlations  attempt 
to  include  anisotropy  and  polarity  by  introducing  a para- 
meter in  addition  to  the  critical  temperature,  pressure  and 
volumes.  The  most  widely  used  is  Pitzer ' s acentric  factor 
(Pitzer  et  al.,  1955,  1957;  Prausnitz,  1969).  other  work 
utilizing  corresponding  states  ideas  includes  fitting 
thermodynamic  properties  (phrased  in  terms  of  integrals  over 
direct  correlation  functions)  to  polynomials  in  density, 
temperature,  and  composition  while  exploring  the  variation 
in  the  coefficients  from  a corresponding  states  viewpoint 


(Mathias,  1978i  Campanella,  1984;  Huang,  1986;  Wooley, 

increases,  the  classical  corresponding  states  correlations 
require  additional  reduced  parameters  (Seed  and  Gubbins, 
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Group  Contribution  Methods 


present  in  molecular  orbital  calculations  (Clark,  1985), 
moments  (Bader  et  al.,  1987;  Bader,  1985),  and  in  integral 


discussion  here  is  limited  to  two  general  classes  of 
semiempirical  group  contribution  methods. 


The  first  approach  supposes  molecular  properties  can  be 
written  explicitly  in  terms  of  the  thermodynamic  state 
variables  and  a set  of  parameters.  It  is  further  assumed 


Doraiswamy,  1965)  and  critical  properties  (Lydersen,  1955). 
The  methods  are  successful  for  compounds  with  rela 


itively 


sadily 


The  second  approach  is  to  assume  that  groups  possess 
thermodynamic  properties  and  molecular  properties  are  given 
as  summations  over  group  properties  (Pierotti  et  al. , 1959). 
The  solution-of-groups  methods  are  an  important  subset  of 
these  methods.  Solution-of-groups  methods  include  the 
analytical  solution-of-groups,  or  "ASOS"  (Wilson  and  Deal, 
1962j  Derr  and  Deal,  1969),  and  UNIFAC  (Fredenslund  et  al., 
1977 ! Macedo  et  al.,  1983).  The  methods  assume  that  group 
contributions  may  be  resolved  into  energetic  and  entropic 
terms,  each  of  which  is  described  by  an  analytic  model. 

engineering  design  calculations  due  to  their  simple  analytic 

Solution-of-groups  methods  such  as  UNIFAC  initially  had 

subsequent  empiricism  has  obscured  the  precise  nature  of  the 
molecular  (and  group)  level  assumptions  at  play.  Attempts 

temperature  dependence  have  proved  difficult  and  have  led  at 

1985).  To  improve  upon  the  approach  in  a systematic  way,  a 
careful  analysis  that  exposes  inherent  assumptions  and 


does  seem  necessary  for  effective  implementation. 


ally  tractable  but  often 


tion.  This  often  limits  use 
redictive  calculations. 


agion  of 
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different  methods  may  lead  to  internal  inconsistencies. 
Group  contribution  methods,  including  the  solution-of- 
groups,  alleviate  these  restrictions  to  some  degree,  but 
certainly  they  can  be  improved  upon.  To  do  so  effectively 
requires  careful  analysis  of  the  methods  followed  by 
modifications,  particularly  those  having  a basis  in 
molecular  theory  that  begin  to  relax  the  more  restrictive 


CHAPTER  III 

STATISTICAL  MECHANICAL  ANALYSIS  OF 
THE  SOLUTION-OF-GROUPS  METHODS 


A basic  assumption  of  solution-of-groups  methods  is 


binatorial  and  residual  contributions.  The  former  includes 


in  solution,  the  residual  chemical  potential  is  written 
within  the  solution-of-groups  method  as 


U^IT) 


- lim  u|$IT)  = 

V*1 


(3-2a) 


and  where  AB  is  the  residual  Helmholtz  free  energy,  NM  is 

types  of  groups  and  C2H5OH  and  C3H7OH  are  different  types  of 
molecules,  though  they  have  the  same  two  types  of  chemically 


chemically  equivalent  groups  are  assumed  to  be  thermo- 
dynamically equivalent ) . 


generality  of  their  conclusions  is  restricted  by  the  well- 

not  apparent  that  all  of  the  molecular-level  physical 
assumptions  required  to  satisfy  eqn.  (3-1)  have  been  exposed 


techniques  based  on  pair  correlation  functions,  where  fully 
rigorous  expressions  can  be  explored,  especially  those  which 
insure  that  the  operation  of  eqn.  (3-2b)  is  meaningful.  In 
particular,  the  analysis  takes  into  account  the  influence  of 
nonnearest  neighbors  and  avoids  the  difficulties  of  lattice 
statistics . 

Here,  the  definition  of  groups,  either  as  groups  of 
atoms  or  some  arbitrarily  defined  interaction  sites,  assumes 
that  the  intermolecular  potential  energy  is  the  sum  of  site- 
site  terms  with  the  site-site  potential  depending  only  on 
the  intersite  distance.  Therefore,  use  is  made  of  the  terms 
"groups"  and  "sites"  interchangeably.  While  the  phrase 
"solution  of  sites"  is  used,  and  eqns.  (3-1)  and  (3-2b) 
imply  that  the  sites  have  properties,  this  should  not  be 
interpreted  to  mean  that  each  site  can  move  independently 
throughout  the  whole  configuration  space.  Such  a picture 
would  mean  that  the  system  would  behave  like  a mixture  of  ns 
species  (where  ns  is  the  number  of  sites)  rather  than  of  nH 
components  (Perry  et  al.,  1981).  This  would  violate  the 
ideal  gas  law  at  low  densities  and  even  allow  for  such 
nonphysical  situations  as  immiscibility,  where  more  sites  of 
a given  type  are  in  a phase  than  the  molecular  stoichiometry 
would  require.  Thus,  the  system  should  be  viewed  as  a 
collection  of  "physical  clusters"  (Chandler  and  Pratt,  1976) 
allowing  for  independent  site-site  interactions,  yet 


constrained  in  configuration  space  so  that  the 
components  (molecules)  may  still  be  identified. 

In  developing  the  necessary  statistical  mechanics,  the 
location  of  the  sites  is  given  either  by  a set  of  vectors 
from  an  origin  to  the  individual  sites  or  by  a set  con- 
taining vectors  from  the  origin  to  the  molecular  centers  of 
mass,  plus  those  from  the  centers  of  mass  to  the  sites  on 
the  molecules.  The  two  sets  span  the  same  configuration 
space  and  both  will  be  used  below. 

There  are  at  least  two  ways  to  perform  statistical 
mechanical  ensemble  averaging  with  standard  classical 
statistical  mechanical  assumptions  (Gray  and  Gubbins,  1984). 
One  could  integrate  over  the  site  position  vectors,  which 
involves  constraints  accounting  for  the  fact  that  the  sites 
are  clustered  together  to  form  molecules.  These  constraints 
are  difficult  to  formulate.  The  other  averaging  procedure 
uses  integrations  over  molecular  orientations  and  center  of 
mass  positions.  This  is  preferable  here,  even  though  a 
formalism  that  is  explicit  in  the  sites  is  sought,  since 
relationships  between  the  site  position  and  the  molecular 
center  of  mass  and  orientation  are  easily  written.  A 
technique  (Ladanyi  and  Chandler,  1975)  that  allows  for 
transformation  of  integration  over  molecular  coordinates 
into  integrations  over  site  positions  is  utilized. 

The  statistical  mechanics  in  terms  of  site-site 
quantities  is  developed  first,  culminating  in  an  expression 


where  <U>  denotes  the  ensemble  averaged  configurational 
energy,  the  quantity  of  interest.  Below  and  in  Appendices  A 
and  B,  derivations  are  given  of  the  rigorous  expression 
relating  <U>  to  interactions  between  sites  of  different 
types  on  the  same  (intra)  and  different  (inter)  molecules. 
Approximations  which  lead  to  solution-of-groups  equations 

Consider  an  arbitrary  space-fixed  coordinate  system 
where  the  total  inter-  and  intramolecular  potential  energy 
functions  are  of  the  site-site  form  and  rewrite  the  ensemble 
averages  over  molecular  center  of  mass  and  orientations  in 
terms  of  integrals  over  the  vectors  from  the  origin  to  the 
individual  sites.  Then,  consider  the  intramolecular  and 
intermolecular  contributions  to  <U>  separately.  These 
techniques  have  been  used  (Chandler,  1982)  to  develop  a 
perturbation  theory  for  the  Helmholtz  function,  but  not  to 
explicitly  treat  the  configurational  energy  as  is  done  here. 

For  any  system  configuration,  the  intermolecular 
potential  energy  can  be  written  schematically  as 


(all  pairs  of  R 
sites  on  different 
molecules ) 


where  u is  the  site-site  pair  potential  function,  and  raw  is 
the  vector  from  the  origin  to  site  number  o (regardless  of 
type ) on  molecule  w ( regardless  of  type ) . This  can  be 
written  more  explicitly  as 


(3-5) 


In  eqn.  13-5)  the  first  two  summations  are  over  pairs  of 
molecules  of  types  M and  M' ; in  a pure  component  there  is 
only  one  type,  but  mixtures  have  more  than  one  type.  Also, 
the  indices  i and  j run  from  1 to  the  number  of  molecules  of 
type  M and  M' , respectively,  and  the  factor  containing  the 
Kronecker  deltas  ensures  that  the  sums  are  over  different 
molecules  (i.e.,  excludes  intramolecular  contributions). 

The  last  sums  are  over  all  sites  (regardless  of  type) 
running  from  1 to  the  number  of  sites  on  a given  type  of 
molecule  (e.g.,  vM).  Finally,  raiM  denotes  the  vector  from 
the  system  origin  to  site  a on  the  ith  molecule  of  type  M. 

We  now  perform  the  ensemble  average,  which  involves  inte- 
grals over  the  positions  of  the  center  of  mass  and  external 
orientations  of  the  molecules,  giving 


J 0 uaMBM'  (raiM,r|3jM' 1 5 


Introducing  the  Dirac  delta  function,  6,  this  is  formally 


Removing  all  nondynamical  variables  from  the  ensemble 
average  yields 


Since  the  molecules  of  a particular  type  (e.g.,  type  M)  are 
indistinguishable,  the  general  relationship 

\l1  6<r  - W 5 = - W> 

is  used  in  order  to  write 


!<»„■  - 6p 


,)6(r' 


Here  note  that  - 6^1 ) terms  survive  from  the 
summations.  The  site-site  pair  correlation  function, 
9oM0M'  < is  defined  as 


where  denotes  the  number  density  of  molecules  of  type  k. 
The  quantity  Sombm'  (r,r' ) is  interpreted  as  the  Joint 
probability  density  of  finding  the  0th  site  (type 
unspecified)  on  any  molecule  of  type  M'  at  location  r'  and 
the  ath  site  (type  unspecified)  of  a molecule  of  type  M at 
position  r.  The  average  intermolecular  configurational 
energy  becomes 

<(Jinter  _ 1 _ f.j  f M M'  -*  -* 

” ^ J o 0 aM3M,lr,r  )PMpM'gaM0M' ,c'r  1 

Selecting  the  origin  to  be  at  r'  and  performing  one  of  the 


integrations  yields 


J uaMBH'(t)9aMBM'(rlr2  **  <3*8> 

0 

The  development  of  the  mixture  intramolecular  contribution 
and  pure  component  cases  proceeds  in  a similar  fashion 
{Appendix  A).  An  alternative  derivation . of  eqn.  (3-8)  uses 
9cMBM'*r'r'>  as  a projection  (Cummings  and  stell,  1982)  of 
the  angle  dependent  molecule-molecule  pair  correlation 
function  g^i (Rm-Rm1 , nM , • ) . This  derivation  is  given  in 
Appendix  B. 

Equations  (3-8)  and  ( A— 3 ) are  site-site  equations  to 
which  the  approximations  are  made  to  yield  the  solution-of- 
groups  expressions  for  the  Helmholtz  energy.  What  follows 
makes  explicit  this  connection. 

The  Gibbs-Helmholtz  equation  integrated  to  obtain  the 
free  energy  terms  of  interest 


Treatments  of  this  integration  may  vary.  For  instance,  one 
may  allow  T<)  to  approach  infinity  and  utilize  an  athermal 
mixture  model  (e.g.,  that  due  to  Guggenheim,  1952)  for  the 
A(Tq)/Tq  term,  thereby  introducing  a "combinatorial" 


expression  (Maurer  and  Prausnitz,  1978).  It  is  also 
possible  to  view  the  A®  term  as  a rigid-body  or  repulsive 
effect  given  an  appropriate  temperature  dependence  of  <U> . 
In  either  case,  collecting  the  inter-  and  intramolecular 
contributions  to  U from  equations  (3-8)  and  (A-3)  gives 

AR1T)  ar(t0) 


Vh- 


where  s^gjjlr)  denotes  the  intramolecular  site-site  pair 
correlation  function  for  molecule  type  M as  defined  in 
Appendix  A. 


The  Solution-of-Groups  Concept 

To  satisfy  eqns.  (3-1)  and  (3-2),  one  must  express  the 
free  energy  in  terms  of  chemically  equivalent  site  types. 
This  is  required  in  order  for  the  operation  of  eqn.  (3 -2b) 
to  be  meaningful.  Presented  here  are  the  most  straight- 
forward yet  plausible  relations  which  accomplish  this, 
namely  assuming  that  the  individual  terms  in  the  summation 


of  component  pairs  in  eqn.  (3-9)  can  be  expressed  in 
specific  ways.  For  the  intermolecular  energy  terms,  for  a 
pair  M,  M',  the  approximations  are 


' I uABtr,gAMBM' (r,r2d 


In  these  equations  the  summations  over  A and  B are  sums  over 
sites  of  chemical  type  A and  B,  respectively,  and  g^lr)  is 
the  intermolecular  site-site  pair  correlation  function  with 
the  site-type  dependence  appearing  explicitly.  There  are 
two  ways  to  deal  with  each  term  in  the  summation  over  M in 
the  intramolecular  terms  of  eqn.  (3-9).  One  is  to  assume 
that  the  intramolecular  contributions  in  the  mixture  are  the 
same  as  in  the  pure  component  case  ( SoMBM  * soMBM>  and  then 
make  site-type  approximations 


ii  “< 


aMBM  'aMBH 


T) r*dr 


(3-lla) 


Without  a composition  dependence  of  s^^lr),  the  AR(To) 
need  have  no  intramolecular  contributions  as  they  would 
always  be  eliminated  in  eqn.  (3-1).  The  alternative,  and 
more  usual  treatment,  is  to  assume  that  SambmIf)  is  indepen- 
dent of  temperature,  so  the  second  integral  of  eqn.  (3-9)  is 
zero,  and  the  intramolecular  composition  dependence  is  put 
in  the  combinatorial  term  of  AR(T0),vlz, 


” J £ VAMVBM  / UAB<rlsAMBMlr,r2  dr  (3-llc) 

Both  of  these  satisfy  eq.  (3-1),  as  is  shown  in  Appendix  C. 

It  is  possible  to  articulate  several  physical  assump- 
tions implicit  in  eqns.  (3-10)  and  (3-11).  Equations 
(3-10a)  and  (3-llc)  can  be  obtained  if  two  assumptions  are 
made.  First,  the  potential  between  all  A and  B type  sites, 
UAMBM ' l r ) * is  written  in  the  form  uM ( r ) , explicit  only  in 


models  of  liquids  treated  intermolecular 


interactions  in  a geometrically-averaged  way  (Barker,  1964). 
An  alternative  physical  assumption  is  that  there  is  an 
effective  set  of  site  "locations”  which  show  sufficient 
symmetry  to  yield  the  correct  total  contributions  to  the 
energy  when  the  pair  parameters  are  obtained  by  fitting 
data.  Molecular  dynamics  simulations  (Lee  et  al.,  1986, 
described  by  Haile,  1986),  on  properties  and  structure  of 
systems  with  one  and  two  site  molecules  seem  to  indicate 
that  some  smearing  may  be  reasonable.  Besides  the  above 
assumptions,  the  choice  between  eqns.  ( 3— lib)  and  (3-llc) 
must  be  made  consistent  with  assumed  temperature  and 
composition  effects  on  the  conformational  changes  upon 
mixing . 

Finally,  there  are  subtle  but  major  effects  of  the 
required  step  of  dropping  the  molecular  labels  (M,  M' ) on 
the  intermolecular  pair  correlation  function  between  eqns. 
(3-10a)  and  (3-10b).  If  the  molecular  labels  were  retained, 
the  use  of  eqn.  (3-3)  is  precluded  and  the  operation  of  eqn. 
(3-2b)  yields  zero.  The  procedure  of  Appendix  C shows  that 
if  the  labels  are  retained,  eqn.  (3-12)  is  required  in  order 
to  satisfy  eqn.  (3-1): 


1/Tq 


1-13) 


uM  ( r ) ( 


where  the  characteristic  length  parameter,  Ojg,  associated 
with  a site  pair  is  apparently  independent  of  the  type  of 
host  molecule,  M'  or  M",  and  the  density  scaling,  o^,  is 
independent  of  the  site  pairs  involved  and  shows  no  explicit 
molecule  pair  dependence.  Note  that  in  arriving  at  eqn.  (3- 
13),  the  dropping  of  molecular  labels  on  the  pair  potential 
involves  a different  assumption  than  suppressing  them  on  the 
correlation  function.  If  eqn.  (3-13)  is  accurate,  the 
intermolecular  contribution  of  eqn.  (3-9)  is 


which  will  satisfy  eqn.  (3-1).  The  scaling  reflected  in 
eqns.  (3-13)  and  (3-14)  is  akin  to  the  n-fluid  theories 
(Leland  et  al. , 1968;  McQuarrie,  1976)  previously  proposed 
for  simple  fluids  and  the  mean  density  approximation 
(Mansoori  and  Leland,  1972;  Ely,  1986b).  Simulation  results 
(Lee  et  al.,  1986,  described  by  Haile,  1986)  indicate  that 


■lithin 


explicit 
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where  < ^.n^,  indicates  a normalized  average  over  molecular 

orientations.  Thus  for  each  value  of  the  molecule-molecule 
center  of  mass  separation  distance,  all  molecular  orienta- 
tions are  sampled  with  their  contribution  to  the  ensemble 
average  weighted  by  the  orientation-dependent  pair  correla- 
tion function.  If  the  molecular  potential  is  of  the  site- 
site  form,  then  eqn.  (4-1)  becomes 

<U>  = V3  M M'  I**"  i®*1'  b uaM0M ^ raM ' r0M ' * ' 

3H,M'«BMA-,aM-QM'»>nB,nM1  (4-2) 

In  the  process  of  performing  an  orientational  average  about 
the  center  of  mass,  the  individual  sites  physically  cover 
the  surfaces  of  spheres.  Thermodynamically  equivalent 
sites,  i.e.  chemically  equivalent  sites  interchangeable  by 
point-group  symmetry  operations,  must  be  the  same  distance 
from  the  molecular  center  so  that  more  than  one  site  may 
sweep-out  or  cover  the  same  spherical  shell.  The  essence  of 


inside 


component  having  a different  shell  radius,  all  were  limited 
to  one  shell  per  molecule.  The  extension  of  the  spherical 
shell  idea  to  multiple  shells  has  been  suggested  (Reed  and 
Gubbins,  1973),  but  has  apparently  not  been  adequately 
developed,  especially  within  a group  contribution  context. 
Before  discussing  the  details  of  uniformly  "smearing" 

describing  site-site  interactions  will  be  examined. 


The  definition  of  the  interaction  sites  is  somewhat 
arbitrary,  although  choosing  them  to  be  the  common  organic 
functional  groups  often  leads  to  accurate  predictions. 

Finer  distinction  between  the  sites  ultimately  results  in 
each  atom  constituting  an  interaction  site.  While  this  may 


reduce  the  total  number  of  types  of  sites,  it  greatly 


increases  the  number  of  interactions  involved  in  a real 
system.  In  addition,  the  origin  of  intermolecular  forces  is 


would  have  to  be  formulated  separately.  At  the  other  limit, 
individual  molecules  could  be  the  interaction  sites,  losing 


tees  typically 


for  ir 


between  small  molecules.  The  model  parameters  become  those 
for  sites  rather  than  the  usual  molecular-level  parameters. 
To  describe  the  potential  energy  between  nonpolar  sites  A 
and  B,  Ujj(r),  spherical  symmetry  is  assumed  and  the  Mie  n-6 
form  (Reed  and  Gubbins,  1973)  is  chosen: 


where  r is  the  site  center-to-center  distance,  oAB  is  a 
characteristic  site  size  parameter,  eAB  reflects  the  depth 
of  the  potential  well,  and  n is  an  integer  (typically, 

9 S n < 30).  The  size  and  energy  parameters  for  the  unlike 
pair  A and  B can  be  related  to  like-like  interaction  para- 
meters by  use  of  mixing  rules,  the  most  common  of  which  are: 

°AB  = l°AA  + °bb>/2.  (4-Sa) 

eAB  = <sAAeBB>*  ( 4-5b) 


Other  model  potentials  for  nonpolar  interactions  (exp-6  and 
Kihara)  are  available  but  are  not  utilized  in  this  work. 

The  form  chosen  correctly  describes  the  r-6  decay  of  the 
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UABlr>  = uM(r,+uAB(r,+uM(r,+uM(r,+uM(c,+uM,r)  (4-6> 


the  u£g(r)  and  ujg(r)  terms  are  zero. 


independent  of 

distribution  factors  into  individual  distributions.  Since 
the  distributions  are  uniform  on  their  respective  shells, 
these  probability  densities  are  constants,  f,  which  can  be 
determined  from  the  normalization  criterion.  For  shell  o or 


Similarly  for 


Bonn': 


probability  density  would  b 


distributed 


shells  the  probability  density  is  again  a constant  and  egns. 
(4-8)  and  (4-9)  simply  become 


where  is  the  number  of  sites  on  shell  a of  molecule  M. 
Note  that  here  the  normalization  criteria  are  somewhat 
different  since  now  there  is  a certainty  of  finding  a site 
within  a surface  area  element  smaller  than  the  total  surface 
area.  Talcing  this  into  account  by  appropriate  modification 
of  the  limits  in  eqn.  (4-7),  one  obtains  (4-10)  and  (4-11). 
consider  a site  at  point  P on  shell  o of  M interacting 

is  b (Figure  1).  The  distance,  r,  from  P to  the  center  of 
some  differential  surface  element,  dAppj, , on  M'  is  deduced 
from  geometric  considerations: 


z 


(0,0, b) 


Figure  1.  Geometry  of 
and  shell  6 


the  interaction  between  p 


For  a site  type  A at  P and  type  B distributed  over  shell  B. 
the  potential  between  P and  dAg^1  is 


u'(r)  s Uftelrl  fgM.  dAg^'  (4-12) 

where  uAB ( r ) is  generally  a sura  of  terms  of  the  form 

uAB<r)  » C • r"n  (4-13) 


with  C a constant  and  ni6.  Appropriate  constants  are 
obtained  from  table  1 and  eqn.  (4-4).  The  net  interaction 
between  P and  shell  B of  M’  is  obtained  by  substituting 
(4-13)  into  (4-12)  and  integrating  over  the  surface  of  shell 
B while  keeping  b fixed.  Equation  (4-12)  then  becomes 


-J  f 


in*. 


where  use  has  been  made  of  eqns.  (4-9)  and  (4-12).  Now 
define 

t s l0M.+b2-2lgM.b  cos«g 


!BM'b  sin*jj  dtg 


integral 


analytically 


to  give 


The  potential  energy  between  shells  is  obtained  by 

over  those  angles  weighted  by  f^.  Exploiting  the  symmetry 
of  the  problem  one  finds: 

b = (laM+RMl'‘21aMRMM'cos  *0]* 


itroducing 


gives  the  net  potential  as 


uaMBM,1RMM')  = J J u'(r,’faM‘*oMdeasin®ad 

= - J u'(r)  sin*ad*a 


The  integral  in  eqn.  (4-17)  c 


aMtBM,RMM' t2"n) 
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where  the  notation  has  been  condensed  by  defining 
R1  3 RMMI  + *aM  + *BM' 
r2  3 rMM'  + loM  " *BM’ 
r3  5 rMM’  " *aM  + *BM' 
r4  3 SMM'  ‘ *oM  ~ *PM' 

Allowing  for  more  than  one  site  per  shell,  multiple  shells 
per  molecule,  and  considering  all  possible  terms  in  eqn. 
(4-6)  with  the  appropriate  values  of  C,  eqn.  (4-18)  can  be 
generalized  to  give  the  net  molecule-molecule  potential. 
The  molecular  pair  potential  is  then  the  sum  of  the  shell- 


In  egn.  (4-19),  the  sums  are  over  the  number  of  shells  on  M 
and  M' , while  the  subscripts  on  the  molecular  center-to- 
center  distance,  R,  are  dropped  here  and  in  subsequent 
equations  for  notational  convenience.  Precisely  which 
chemical  site  type  occupies  each  shell  is  understood  to  be 
user-supplied.  For  a given  chemical  site  type,  not  all  the 
terms  in  (4-19)  will  be  present,  but  they  are  included  in 
this  and  subsequent  derivations  for  the  sake  of  generality. 
Once  again,  if  net  molecular  moments  are  calculated  from 
site  moments,  the  form  of  the  polar  contributions  to  (4-19) 
remains  the  same  except  that  the  site  a,  u and  Q are 
replaced  with  the  molecular  analogues. 


One  form  of  egn.  (4-19J  useful  for  exploring  properties 
of  the  shell  potential  with  variations  in  the  parameters  is 
that  applicable  to  nonpolar  sites,  that  is. 


which  will  be  used  below.  Before  beginning  a detailed 
discussion  of  the  properties  of  the  shell  potential,  it 
should  be  noted  that  for  the  case  of  one  shell  per  molecule, 
one  site  per  shell  and  equal  shell  radii  (pure  components), 
eqn.  (4-20)  matches  exactly  the  potential  used  in  earlier 
virial  coefficient  calculations  (DeRocco  and  Hoover,  1962). 
For  unequal  shell  radii  [binary  mixtures),  again  a previous 
result  (Dorsch  and  DeRocco,  1971)  is 


For  the  case  of  Lennard-Jones  sites  this  becomes: 


recovered. 


,,-*0  gives 


lint 

PM' 


Again  using  l'Hopital's  rule,  one  finds: 


[(R+lPM 


which  matches  eqn.  (4-13)  exactly. 

Careful  examination  of  the  shell-shell  potential,  e.g. 
eqn.  (4-18),  indicates  that  when  R approaches  , the 


potential  becomes  infinite  due  to  R4  being  vanishingly 
small.  Thus  the  potential  has  a spherical  hard  core 
repulsion  and  thus  is  similar  in  form  to  the  spherical  core 
Kihara  potential  (Kihara,  1976).  For  realistic  inter- 
actions, the  Kihara  potential  assumes  that  the  molecular 
pair  potential  depends  only  on  the  shortest  surface-to- 
surface  distance,  d,  between  the  repulsive  cores.  The 
actual  interaction  is  given  by  the  chosen  pair  potential 
model  evaluated  at  d.  The  shell-shell  model  proposed  above 
(e.g.  eqn.  4-18)  evaluates  the  potential  not  only  between 
the  two  closest  points,  but  at  distances  between  all  points 
on  the  surfaces  of  the  cores  or  shells.  One  might  question 
whether  or  not  the  shell-shell  potential  can  be  rewritten  in 
the  same  functional  form  as  the  Kihara  potential.  Using  a 
potential  of  the  form  of  eqn.  (4-13)  with  equal  core  radii, 
the  Kihara  form  can  be  written  as 

UMM'<R>  = C'-(R-2i)'n  (4-22) 

where  c’  is  a constant.  The  actual  value  of  the  constant  is 
not  important  since  the  essence  of  the  problem  is  to 
determine  if  eqn.  (4-18),  with  equal  shell  radii  and  single 
site  types,  can  be  rearranged  into  the  functional  form  of 

u"m,(b’  = 2^'I^'7I^[(R+2t)3’n'2R3"n+lB'21)3 


!"n](4-23) 


be  recasC  into  the  form  of  eqn.  (4-22)7  Considering  the 
case  of  n=4,  eqn.  (4-23)  becomes: 


C rR(R-2i)-2(R+2i)(R-2i)+R(R+2t)' 
8i2R  V,  (R+2i)R(R-2t) 


n be  simplified  to  read 


Equation  (4-24)  is  not  of  the  form  of  (4-22)  since 


Consequently,  one  must  conclude  that  the  shell-shell 
potential  is  different  than  the  spherical  core  Kihara  model 
for  pair  potentials  of  the  form  considered  in  this  work. 

Nonspherical  repulsive  cores  may  be  used  in  the  Kihara 
potential.  However,  such  cases  require  two  Cartesian 
coordinate  systems  with  parallel  axes  and  movable  origins 
located  within  the  molecules  (Boublik  and  Nezbeda,  1986). 

The  calculation  of  properties  is  then  a non-trivial  task  due 
to  the  fact  that  d is  now  orientation  dependent. 


spherically  symmeti 


the  molecular  Lennard-Jones  potential  is  shown  in 
Figure  2.  The  potentials  were  adjusted  so  as  to  have  a 
common  minimum.  The  figure  indicates  that  while  similar 
in  form,  the  shell  potential  has  a narrower  well.  The 
effect  of  an  additional  shell  on  the  location  and  depth  of 

shells  use  Lennard-Jones  site  interactions.  Depending  on 
the  radius  of  the  inner  shell,  the  net  molecule-molecule 
potential  can  vary  significantly  both  qualitatively  and 
quantitatively  given  a fixed  larger  radius.  These 


form  are  included,  and  when  additional  shells  are  present. 
Thus  the  concentric  shell  model  can  yield  a 
wide  range  of  spherically  symmetric  molecular  potential 
energy  functions.  The  figure  also  indicates  how  aspects 

the  potential  since  for  different  shell  radii,  the  inner 


shell  has  different  effects  on  molecule-molecule  potential. 
The  influence  of  this  shell  varies  with  shell  radius. 


us,  consider  the  central 
er  shell  would  likely  not 


radii 
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simulations  using 


potentials  confirm 


( Lustig, 


A central  idea  of  group  contribution  methods  based  on 

chemical  type  on  a given  host  molecule  are  treated  as 
thermodynamically  distinct  when  they  cannot  be  interchanged 
by  using  the  point  group  symmetry  operations.  Thus,  if  two 


chemically  equivalent  site  types  are  thermodynamically 
equivalent.  The  success  of  these  methods  for  such  compounds 
as  linear  hydrocarbons  indicates  that  the  CH3  and  CH2  groups 


Figures  4a  and  4b  compare  the  shell-shell  potential  for  a 


tions  with  the  potential  obtained  from  a single  shell 
analogue  occupied  by  three  sites  with  identical  Lennard- 


similar. 
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from  eqn.  ( 4-5b>  while  o is  determined  from  (4-5a)  with 
slight  modification: 


This  accounts  for  more  than  one  site  per  shell.  The  strict 
microscopic  comparison  of  potentials  indicates  reasonable 
agreement.  The  differences  can  likely  be  attributed  to 
inadequacies  in  the  mixing  rules.  Mixing  rules  often 
require  empirical  modification  to  yield  accurate  corre- 
lations (Prausnits,  1969).  Thus  it  can  be  concluded  that 
the  shell  model  is  not  generally  incompatible  with  three 
parameter  corresponding  states  correlations  provided  that 
all  site-site  interactions  are  conformal.  Note  that  for 
large  differences  in  a and  e,  the  one-shell  analogue  can  be 
somewhat  inaccurate  for  conformal  potentials  as  shown  in 
Figure  7 . The  discrepancy  here  is  probably  due  to  a 
breakdown  of  three  parameter  corresponding  states  theory 
itself.  Treatment  of  larger  molecules  as  well  as  different 
types  of  site-site  interactions  (e.g.  polar  sites) 
diminishes  the  usefulness  of  a three  parameter  description 
and  four  or  more  parameters  are  then  required  in  corres- 
ponding states  correlations  (Reed 


and  Gubbins,  1973). 
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spherical  shell  model  retaining 


has  been  introduced.  The  model  assumes  that  sites  may  be 
uniformly  distributed  over  the  surfaces  of  spherical  shells 
centered  at  the  molecular  center  of  mass.  For  any  given 
molecule  all  chemically  equivalent  sites  that  are 
thermodynamically  equivalent  will  reside  on  the  same  shell. 
Properties  obtained  from  this  potential  are  dependent  upon 
the  assumption  that: 


I ^ '“m'0!!' lgMM' '*[^,0,^, 

"(j^M  U‘251 

where  the  quantity  on  the  right  hand  side  of  (4-25)  utilizes 
the  shell-shell  potential  and  associated  pair  correlation 
function.  The  shell  radii  are  taken  as  parameters  adjusted 
such  that  eqn.  (4-25)  is  approximately  satisfied.  Pure 
component  data  are  suggested  as  a means  of  deducing  the 

Site-site  interactions  are  assumed  to  be  of  the  Mie  n-6 
form  for  nonpolar  sites.  For  polar  sites  two  schemes  are 
suggested.  The  first  would  utilize  a prescription  for 
adding  site  electrostatic  moments  to  obtain  net  molecular 
moments.  These  would  be  located  at  the  center  of  mass  and 


e-averaged  to  yield  spherically  symmetric 
potentials.  The  second  alternative  would  orientationally 

the  shells.  Again,  spherically  symmetric  interaction 
potentials  are  obtained.  For  both  cases,  the  dispersion  and 


include  nonpolar,  dipolar,  quadrupolar  and  polarizability 
were  developed  for  shell-shell  interactions. 

the  constituent  shell-shell  interactions  was  found  to  have 
well-behaved  properties.  For  single  shell  molecules  with 

recovered  in  the  limit  of  vanishing  shell  radii.  For 
multiple  shell  molecules  occupied  by  chemically  equivalent 
site  types,  the  sites  apparently  can  act  as  if  they  were 
thermodynamically  equivalent  in  specific  cases  and  the 
solution-of-groups  form  can  be  matched.  The  shell  potential 

corresponding  states  correlations  for  "normal"  fluids.  The 
model  can  be  made  quite  flexible  by  introducing  additional 


ideas 


Group  contribution  equations  of  state  (GC-EOS)  have 
been  proposed  to  treat  gas  solubilities  in  liquids  (Skjold- 


ealculate  vapor-liquid  equilibria  (Gupte  et  al.,  1986).  The 
first  of  these  uses  mixing  rules  based  on  the  NRTL  model 


of  empirical  factors  into  the  mixing  rules.  The  third 
state  to  that  of  the  solution-of-groups  method,  UNI FAC,  in 


introduced  to  allow  the  equation  of  state  parameters  to  be 
temperature  dependent.  Careful  examination  of  predictions 


not  correct  (Gupte  et  al.,  1986).  Another  problem  with 
this,  and  related  methods  (Vidal,  1978;  Huron  and  Vidal, 
1979),  is  that  in  the  low  density  limit  these  GC-EOS  yield 


liquid 


parameters  possessing  physical 


where  v = i and  a and  b are 
interpretations.  The  "b"  parameter  represents  the 
incompressible  volume  of  closed-packed  molecules  while  "a" 
is  associated  with  attractive  forces  between  molecules  and 
decreases  the  pressure  below  the  ideal  gas  value.  The 
equation  is  cubic  in  volume,  shows  a vapor-liquid  phase 
transition,  and  a critical  point  (Murrell  and  Boucher, 

1982). 

The  van  der  Waals  form  of  the  equation  of  state  can  be 
derived  from  statistical  mechanical  perturbation  theory. 

Here  such  a derivation  is  given  so  that  the  physical 
significance  of  the  parameters  can  be  expressed  explicitly 
in  terms  of  molecular  quantities.  The  shell-shell  potential 
is  used  in  a generalization  to  mixtures  of  a previously 
given  derivation  (McQuarrie,  1976).  Assume  that  the  system 
potential  energy,  U,  can  be  resolved  into  two  parts: 

U = U<°>  + U<1)  (5-1, 

where  u!°*  is  the  potential  energy  of  a reference  system 
(which  is  designated  below)  and  U^1'  consists  of  the 
perturbation.  Each  term  in  eqn.  (5-1)  is  a function  of  the 
set  of  molecular  center  of  mass  vectors  {R}.  No  angular 
dependence  is  present  due  to  the  form  of  the  shell-shell 
potential.  Furthermore,  the  explicit  dependence  on  (R)  is 
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suppressed  for  notational  convenience.  The  configurational 
integral,  Z,  can  then  be  written  as 


where  Z*°*  is  the  configurational  integral  of  the  reference 
system,  < >g  denotes  a canonical  ensemble  average  in  the 
reference  system,  and  3 = 1/kT.  Assume  that  U(1)  is  small 
relative  to  kT.  The  exponential  term  in  eqn.  (5-2)  can  then 
be  approximated  by 


= Jd<R)exp(-BU(0): 


Jd<R)expt-3(U(0,+UU1n 

|d(R)exp(-pu101) 


<exp(-30*^l )>o  * i-p<u(D >0  + ... 


expl-KPU^lig) 


(5-3) 


which  is  similar  in  spirit  to  the  u-expansion  discussed  in 
Chapter  2.  Equation  (5-2)  now  becomes 

Z = Z<°)exp(-<pU<1)>0)  (5-4) 

When  U<1>  is  the  sum  of  pair  potentials  between  molecules, 
one  can  use  the  energy  equation  to  write: 

<u(1)>0  = ~ E s,  nh  HH'J  uMM!|R)9(m!<!1>R2'iR  <5-S) 

where  g^!(R)  is  the  pair  correlation  function  in  the 
reference  system.  Introduction  of  the  mean-field  approxima- 


gm  (B) 


C 


if  R<RC(MM'  ) 
if  R2SC(HM'  ) 


allows  eqn.  (5-5)  to  be  written  as 


wherein  it  should  be  noted  that  the  parameter  K can  be  a 
function  of  the  state  conditions  and  RCIMM' ) is  a charac- 
teristic cut-off  distance  for  the  pair  M-M' . Choices  for 
these  are  discussed  below.  Now  define  the  quantity 


mole  fraction  of  species 


Substitution  of  egn.  (5-7)  into  (5-4)  gives: 


2 - z(O*exp(0aNo) 


The  thermodynamic  pressure  is  given  by  the  canonical 
relation 


One  could  employ  the  hard  sphere  equation  of  state 
(Carnahan  and  Starling,  1969)  or  the  hard  convex  body 
equation  (Boublik  and  Nezbeda,  1986)  for  P0,  both  of  which 
introduce  a composition  dependence  into  P0 . Here  we  shall 
consider  the  reference  system  of  hard  spheres  whose  poten- 
tial energy  is  infinite  at  contact  and  zero  elsewhere. 


which  upon  substitution  of  eqn.  (5-9)  becomes: 


(5-10) 


An  approximation  to  Z for  such  a system  can  be  obtained  by 
consideration  of  excluded  volume  effects.  Each  molecule  of 
type  M in  the  reference  system  occupies  a volume,  v^,  given 
by: 


where  d^  is  the  diameter  of  sphere  M.  In  a c-component 
reference  system  of  hard  spheres,  the  volume,  VM,  excluded 
to  each  sphere  of  type  M is  approximately: 


where  an  additional  factor  of  j was  introduced  since  for  an 
interacting  pair  M-M' , half  the  excluded  volume  effect  is 
assigned  to  M.  The  reference  configurational  integral,  now 
of  the  form  of  a "free"  volume,  can  be  approximated  as: 


where  V is  the  system  volume.  The  reference  pressure,  Pf 


(5-11) 
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Substituting  eqn.  (5-13)  into  (5-10)  gives 


which  is  the  van  der  Waals  equation  of  state. 

Ideally  one  would  like  to  express  both  the  "a"  and  "b" 
parameters  in  terms  of  group  contributions.  This  has 
explicitly  been  done  for  "a"  through  its  definition  in  eqn. 
(5-8).  Choosing  the  hard  sphere  diameters  in  the  reference 
system  to  correspond  with  the  zero  of  the  shell-shell 
potentials  implicitly  includes  group  contributions  into  the 
"b"  parameter.  Due  to  the  form  of  the  shell-shell  poten- 
tial, this  zero  (which  is  analogous  to  the  Lennard-Jones  o) 
must  in  general  be  located  numerically.  This  choice  for  dfl 
will  be  utilized  in  the  following  chapter.  More  traditional 
choices  for  the  "b"  parameter  are  possible,  such  as  allowing 
it  to  be  an  empirically  adjusted  factor  or  obtaining  it  from 
pure  component  density  data.  In  light  of  the  group 
contribution  formalism  developed  here,  these  options  are  not 
pursued. 

The  choice  of  the  zero  of  the  shell-shell  pair  poten- 
tial as  the  hard-sphere  diameter  allows  RC(MM' ) in 


eqn.  15-8)  to  be  chosen  in  a consistent  fashion.  For  a pair 
M-M'  the  characteristic  cut-off  is  taken  to  be 


(5-16) 


where  d^i  is  the  zero  of  the  pair  potential  between  M and 
M' . Equation  (5-8)  then  becomes 


• - v*  J "Si ,ia*  i"-17' 

Sm1 

Since  the  potential  is  zero  in  the  reference  system  for  all 
> the  quantity  u^J  (R)  is  simply  the  shell-shell 
potential.  The  integrals  resulting  from  the  substitution  of 
eqn.  (4-19)  into  (5-17)  can  be  evaluated  analytically 
(Appendix  D).  The  resulting  expression  for  "a"  is: 


(5-18) 


Sw'kt  5 ldMMl+<~1,k*a 


and  again  note  that  not  all  terras  in  eqn.  (5-18)  need  be 
present  £or  a given  pair  M-M' . 

The  expression  derived  here  for  the  "a"  parameter  is 
quadratic  in  the  mole  fraction  and  as  such  will  yield  the 
correct  composition  dependence  of  the  second  virial  coeffi- 
cient in  the  low-density  limit.  This  limit  is  not  properly 
met  with  the  presently  available  GC-EOS  of  this  form. 

Application  of  the  van  der  Waals  equation  to  mixtures 
in  a non-group  contribution  format  typically  involves  mixing 
rules  for  the  "a"  parameter.  The  most  common  of  these 
(Rasmussen.  1987)  are: 


aMM'  “ (5"‘ 

where  k^i  is  an  empirical  correction.  The  use  of  eqn. 
(5-18)  avoids  the  use  of  eqns.  (5-20)  or  other  empirical 


mixing  rules.  Using  (5-18)  gives  the  a^.  contributions  to 
"a"  directly  as  integrals  over  the  shell-shell  pair 
potential  between  M and  M' . 

The  mean-field  approximation  to  the  pair  correlation 
function  introduced  in  eqn.  (5-6)  can  be  dealt  with  several 
ways.  The  parameter  K is  in  general  a function  of  state 
conditions.  Different  values  of  K for  each  distinct 
molecular  pair  M-M'  could  be  used,  thereby  introducing 
composition  dependence  into  the  mean  field  approximation. 
This  would  substantially  increase  the  number  of  parameters 
involved  and  thus  is  not  examined  in  this  work.  Density 
effects  on  K can  be  introduced  by  consideration  of  density 
expansions  of  the  pair  correlation  function.  This  option 
will  be  discussed  in  the  next  section.  For  a GC-EOS  of  the 
van  der  Waals  form,  the  available  methods  indicate  that 
incorporating  only  a temperature  dependence  into  the  K 
parameter  is  sufficient  for  reliable  predictions.  Thus,  the 
form  of  the  temperature  dependence  of  K will  be  explored 
first. 

Recently  proposed  GC-EOS  (Gupte  et  al.,  1986)  have 
taken  the  "a"  parameter  to  be  temperature  dependent  by 
assuming  the  ad  hoc  form: 

a( T ) - [ 1+Ci ( 1-Tri ) +c2 ( 1-Tr 1 ) 2+C3 ( 1-Tr * ) 3 ] 2 (5-21) 

where  Tr  is  the  temperature  reduced  by  its  critical  value. 


This  fc 


riginally 


more  traditional 


J...J  [X-u/kT)+}(U/kT)2 ]dR3-  -dR,, 

J---J  [1-  /kT+itU/kT)2 IdRj^-'  dRj, 


where  the  factors  and  contain  the  integrals  of  eqn . 
(5-23),  and  assume  the  role  of  the  constants  in  eqn.  (5-21). 
Upon  truncation  of  the  power  series  in  eqn.  (5-24),  the 
resulting  expression  is  of  the  Pade  form.  Pade  approximants 
are  known  to  be  a powerful  means  of  approximating  truncated 
power  series  (Ralston,  1965).  Recently,  ratios  of  poly- 
nomials have  been  proposed  without  justification  to  describe 
the  temperature  dependence  of  the  "a"  parameter  (Dimitrelis 
and  Prausnitz,  1986).  If  a three  parameter  model  for 


temperature  dependent  "a"  parameters  is  required,  an 


is  suggested  for  use  in  eqn.  (5-18). 

Equations  (5-15),  (5-18),  and  (5-24)  along  with 
knowledge  of  the  "b"  parameter,  which  is  directly  related  to 


the  zero  of  the  shell-shell  pair  potentials  through  eqn. 
(5-14),  constitute  a GC-EOS  of  the  van  der  Waals  form.  This 


GC-EOS  can  be  used  to  develop  group  contribution  correla- 
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<n-3)  J 


where  dj^t  has  been  used  as  the  characteristic  cut-off 
distance  so  that  the  notation  of  egns.  (5-19)  may  be 
employed.  The  quantity  "a"  in  eqn.  (5-27)  is  given  by 
(5-18). 

If  K is  taken  to  be  independent  of  density,  the  GC-EOS 
of  eqn.  (5-27)  is  quadratic  in  density  and  as  such  will  not 
show  a critical  point  or  a vapor-liquid  phase  transition. 

To  alleviate  this,  empirically  introducing  density  dependent 
K parameters  or  density  expansions  of  the  pair  correlation 
function  might  be  considered.  Consider  the  latter 
alternative.  As  in  the  discussion  of  the  high-temperature 
expansion  of  the  pair  correlation  function  given  above,  the 
focus  is  on  the  pure  component  case  only  for  notational  and 


algebraic  convenience.  The  conclusions  of  the  discussion 
again  apply  to  the  multicomponent  fluid.  Consider  the 
expansion  of  9mm(R/P«T)  : 

gMMlB.O/T)  = g0(R,T)  + pg^R.T)  + p2g2(R,T)  + •••  15-28) 


where  in  the  mixture  case,  the  g^  are  also  functions  of 
composition.  One  might  then  truncate  eqn.  (5-28)  and 
substitute  into  (5-26)  to  introduce  additional  density 
dependence.  However,  due  to  the  relationship  (McQuarrie, 
1976)  between  the  jth  virial  coefficient,  Bj,  and  the  terms 
in  (5-28): 


Bj+2 


dWR) 

dR 


gjfr.T) 


R3dR 


(5-29) 


substitution  of  (5-28)  into  (5-26)  simply  reproduces  the 
virial  equation  of  state.  Virial  coefficients  for  the 
shell-shell  potential  (5-26)  could  be  obtained  more  directly 
by  use  of  the  formulae  relating  them  to  Mayer  f-functions 
(Reed  and  Gubbins,  1973).  If  empirical  K parameters 
dependent  on  T and  p are  introduced,  it  must  be  done  with 
care  so  that  the  proper  mole  fraction  dependence  for  the 
term  cubic  in  density  is  recovered.  It  should  be  cubic  in 
mole  fraction  so  as  to  match  that  of  the  third  virial 
coefficient  (Hirschfelder  et  al.,  1954).  However,  models 
cubic  in  mole  fraction  do  not  accurately  represent  liquid 


functions.  The  resulting  equation  of 


tility 


emphasized 


CHAPTER  VI 

CALCULATION  OF  BULK  PROPERTIES  FROM  THE 
CONCENTRIC-SHELL  POTENTIAL  ENERGY  MODEL 

The  concentric  shell  model  proposed  and  developed  in 
the  preceding  chapters  ignores  certain  details  of  the 

nonsphericity  and  the  presence  of  preferred  orientations 
between  molecules  are  not  treated  explicitly.  These  effects 
are  incorporated  implicitly  through  the  parameterization 
(e.g.  adjustment  of  shell  radii).  The  intent  of  this 
chapter  is  to  establish  the  feasibility  of  using  the 
concentric  shell  potential  to  calculate  thermodynamic 
properties  of  molecular  fluids.  Both  the  vapor  and  liquid 

An  essential  requirement  for  use  in  a group  contribu- 
tion format  is  the  transferability  of  site-site  potential 

given  chemical  site-type  must  be  independent  of  the  host 
molecule  type.  Transferability  of  site  parameters  as  well 
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the  number  of  chemically  distinct  site-types  present  in  the 
system.  The  group  contribution  formalism  proposed  in  this 
work  includes  at  most  three  site  parameters  per  chemically 


and  an  electrostatic  property  such  as  the  polarizability, 
dipole  moment  or  quadrupole  moment.  The  UNI? AC  solution-of- 


chemically  different  sites.  A comparison  of  the  number  of 
site  parameters  required  in  the  proposed  method  with  the 
number  required  by  UNIFAC  is  given  in  Table  2.  As  the 
number  of  site-types  present  increases,  the  shell-shell 
model  generally  requires  fewer  parameters  than  UNIFAC. 
Furthermore,  if  the  electrostatic  effects  are  ignored  or 

type.  Determination  of  site  size  parameters  (i.e.  o)  from 
the  number  of  adjustable  parameters  required  by  the  shell 


o the  number  require 


The  globular  molecule  0(01)3)4  was  considered  first 
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COMPARISON  OF  SITE  PARAMETER  FOR 
THE  UNIFAC  AND  SHELL  MODELS 


of  chemically 


6 2 

12  12 

IS  20 


(5-15),  < S— IS ) and  (5-18).  All  site-site  interactions  were 
assumed  to  be  of  the  Lennard-Jones  form.  The  central  carbon 
atom  was  placed  on  a shell  of  radius  0.00001A.  The  four  CH3 
groups  were  distributed  over  a shell  of  radius  l.SA,  which 
is  approximately  the  carbon-carbon  bond  length.  The  o and  e 
values  for  the  083  group  were  taken  to  be  adjustable.  The  a 
values  for  the  CH3  and  central  C groups  were  chosen  such 
that  their  ratio  was  0.67  which  roughly  matches  the  ratio  of 
group  diameters  as  deduced  from  van  der  Waals  group  volumes 
(Bondi,  1968).  In  a similar  vein,  the  ratio  of  e for  the 
CH3  group  to  that  for  the  central  C group  was  chosen  to  be 
0.37.  Specifying  these  ratios  a priori  served  to  reduce  the 
dimension  of  the  parameter  space  in  which  the  search  for 
optimal  parameter  values  took  place.  The  o and  e parameters 
for  the  interactions  between  the  CH3  and  C sites  were 
obtained  with  the  mixing  rules  eqns.  (4-5).  Furthermore, 
the  mean  field  parameter,  X(T),  which  approximates  the  pair 
correlation  function  was  taken  to  be  adjustable. 

Table  3 contains  calculations  for  the  pressure  and 
Helmholtz  energy  departure  functions  for  C(CH3)4.  These  are 
compared  to  the  experimental  PVT  data  (Das  et  al.,  1977). 

The  Helmholtz  departure  function,  AD,  is  defined  to  be  the 
Helmholtz  energy  minus  the  Helmholtz  energy  of  an  ideal  gas 
at  the  same  temperature  but  at  the  standard  state  of  one 
atmosphere  pressure.  The  experimental  values  of  the  volume 
and  temperature  were  utilized  in  the  GC-EOS  to  calculate  the 


when  using  parameters  deduced  £rom  a different  region  of  the 
phase  diagram  (Abbott,  1979). 

same  values  for  the  parameters  o and  e to  determine  if  the 
saturation  conditions  were  accurately  reproduced.  The  mean 
field  parameter  value,  K = 2.27,  accurately  reproduced  the 
vapor  phase  pressures  (less  than  1%  error)  and  free  energies 
(less  than  1%  error)  up  to  the  saturated  vapor  conditions. 
Using  the  saturated  vapor  volume  as  input  gave  a calculated 
pressure  of  10.26  atm.,  versus  the  experimental  value  of 
10.36  atm.,  while  at  the  saturated  liquid  condition  the 

value  of  10.26  atm. 


The  inability  tc 
envelope  is  likely  dt 


duce  vapor-liquid  pha 


adequately  repro- 


aanssaad  aaoncna 


to  give  proper  vapor-liquid  phase  behavior,  the  equations 
often  do  not  accurately  reproduce  liquid  phase  densities. 

The  fact  the  e/k  is  so  small  (giving  essentially  hard  sphere 
interactions)  suggests  that  the  parameter  values  are  not 
optimal.  Predictions  from  the  GC-EOS  proved  quite  sensitive 
to  the  values  of  the  parameters.  For  example,  using  the 
above  values  which  reproduce  the  saturation  conditions  along 
the  T = 430°K  isotherm,  a 1%  increase  in  e/k  gives  rise  to  a 
10.6%  decrease  in  the  calculated  pressure  at  the  saturated 
liquid  volume.  A 1%  increase  in  the  value  of  the  CH3  a 
value  leads  to  a 6.4%  decrease  in  that  calculated  pressure, 
while  a 1%  increase  in  the  value  of  K results  in  a 10.6% 
decrease  in  the  calculated  pressure.  Much  smaller  effects 
on  the  calculated  free  energies  and  saturated  vapor 
pressures  were  observed.  This  sensitivity  to  parameter 
values  prevented  accurate  representation  of  the  entire  phase 
envelope. 

The  calculations  presented  here  for  0(013)4  indicate 
that  the  GC-EOS  can  reproduce  gas  phase  data  accurately  and 
can  be  forced  to  match  specific  saturation  data  by  careful 
adjustment  of  the  parameters.  However,  optimization  would 
require  a larger  database  along  with  careful  selection  of 
initial  parameter  guesses.  Since  the  purpose  here  is  to 
establish  credibility  for  the  proposed  method,  such 
optimization  has  not  been  done.  Further  calculations  will 
be  required  to  explore  the  full  transferability  of  the  site 
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Table  4 

PRESSURES  CALCULATED  FROM  THE  VAN  DER  WAALS 
GC-EOS  FOR  PROPANE 

T = 560°K  IK(T)  = 0.001) 


1 (KIT)  = 0.00128) 


0.9874 

4.911 


14.93 

24.62 


-5. 24X10-2 
-4.78X10"1 
5.50x10“* 


-2.27X10"1 


PRESSURES  CALCULATED  FROM  THE  VAN  DER  WAALS 
GC-EOS  FOR  BUTANE 


T = 510.93°K  (KIT)  = 0.58) 
P-calc.  (atm)  P-exp.  (atm) 


T = 477°K  (K(T)  = 0,615) 
P-exp.  (atm) 


9.869 


29.608 


1.94x10-1 


-3. 24X10'1 
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continually  rescaled  to  maintain  the  specified  initial 
value.  When  the  mean-squared  displacement  of  the  particles 
reached  a value  of  0.4o,  the  temperature  scaling  ceased  and 
the  equilibrium  portion  of  the  simulation  began.  The  Verlet 
neighbor  list  (Verlet,  1967)  was  used  as  a time  saving 
device  as  was  the  minimum  image  criterion  (Wood,  1968).  The 
shifted  force  potential  (Nicolas  et  al.,  1979)  was  used  with 
a cut-off  distance  of  2. So.  Long-ranged  corrections  to  the 
shifted  force  results  (Nicolas  et  al.,  1979)  were 
incorporated.  Details  of  these  techniques  as  well  as 
calculation  of  bulk  properties  from  the  phase  space 
trajectories  are  contained  in  the  literature  cited  and  in 
textbooks  on  the  subject  (Heermann,  1986).  The  number  of 

to  run  in  this  study,  thus  the  number  of  steps  will  be 
reported  along  with  the  simulation  results.  The  program  was 
initially  tested  by  matching  results  of  Lennard-Jones  fluids 
by  simulating  concentric  shell  fluids  with  a shell  radius 
approaching  zero.  The  timestep  was  always  adjusted  to 
insure  that  the  total  energy  was  conserved  to  five  places 
throughout  the  course  of  the  run.  The  kinetic  and  potential 
energy  were  monitored  and  observed  to  oscillate  180°  out  of 
phase,  as  should  be  the  case.  A FORTRAN  program  for 
performing  molecular  dynamics  simulations  for  particles 
interacting  through  eqn.  (4-20)  is  given  in  Appendix  G. 
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center  Lennard-Jones  potential  energy  model  (Cheung  and 

diatomic  simulations,  the  interpolation  formula  proposed  by 
Wojcik  et  al  (1982)  was  utilized  to  yield  comparisons  at  the 
same  temperature.  This  interpolation  formula  reproduces  the 


hermodynamic  ener 


The  results  obtained  from  simulations  of  concentric 

simulations  in  Table  6.  These  are  reported  at  three  reduced 
densities,  p*  = po^ , and  five  reduced  temperatures, 

T*  = kT/e.  The  densities  were  chosen  to  match  those  used  in 
the  nitrogen  simulations  while  the  reported  values  of  T*  are 


orted.  The  agreement  of  ther: 
orable.  The  values  of  the  re 


tently  lower.  The  pressure  proves  to  be  a difficult 
quantity  to  calculate  from  simulation  due  to  the  fact  that 
the  virial  function  is  continually  changing  signs.  This 


occurs  since  at  liquid- like  densities,  molecules  tend 
congregate  about  the  minimum  of  the  potential  energy. 
Cheung  and  Powles  simulated  the  diatomic  nitrogen  for 
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timesteps  and  obtained  good  statistics  for  the  average 
pressure,  but  due  to  the  steepness  of  the  concentric  shell 
potential  about  the  minimum  (see  Figure  2)  the  fluctuations 


runs  may  be  required  to  obtain  reliable  values  for  the 


values  of  the  energy  and  its  variation  with  state  condi- 
tions, all  other  thermodynamic  quantities,  in  principle,  can 


tion,  obtaining  accurate  pressures  from  simulation  was  not 


pursued  in  this  study. 


The  simulation  results  summarized  in  Table  6 indicate 
that  the  concentric  shell  model  can  be  used  to  model  dense 


results  obtained  here  suggest  further  investigation. 
Methods  other  than  simulation  for  investigating  this 


The  use  of  computer  simulation  to  calculate  bulk 
properties  for  dense  fluids  proves  to  be  an  ambitious 
project  even  for  concentric  shell  molecules.  While  some 
degree  of  correlation  as  with  the  dumbbell  fluids  is 


du Id  probably 


be  required  for  each  "cla 

molecules  with  a certain  number  of  shells  (i.e.,  thermo- 
dynamically distinct  sites).  In  addition,  within  each  class 
different  correlation  and  parameterization  would  likely  be 

interactions  such  as  dipole,  quadrupole,  and  nonpolar 
interactions.  In  light  of  this,  as  well  as  the  general 

liquid  region  of  the  phase  diagram,  other  treatments  of 
dense  fluids  based  on  the  concentric  shell  potential  can  be 


The  form  of  the  concentric  shell  potential  is  similar 
to  the  Lennard-Jones  potential  when  the  latter  is  adjusted 
so  that  the  two  well  minima  coincide  (Figure  2).  As  a 
consequence,  the  first  and  major  peaks  of  the  pair  correla- 
tion functions  overlap  to  a significant  degree.  Use  of 
appropriately  scaled  Lennard-Jones  pair  correlation  func- 
tion for  the  pair  correlation  function  is  worthy  of 
consideration.  Here  the  perturbation  parameter  is  the 


pair  potentials.  The  formal  perturbation  expansions  for  the 
pair  correlation  function  are  available  (Gray  and  Gubbins, 
1984).  In  addition,  pair  correlation  functions  for  the 
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shell. 


considered . 
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These  events  are  independent  due  to  the  assumption  of 
uniform  site  distributions.  To  relate  the  two  correlation 
functions,  all  values  of  the  center-to-center  distance  are 
allowed  as  are  all  values  of  the  angles  eA,  eB,  *A  and  $B. 
Each  is  weighted  by  the  probability  of  that  configuration 
and  only  those  configurations  having  shells  touching  the 
specified  points  are  retained.  This  "filtering”  is 
accomplished  once  again  by  the  introduction  of  Dirac  delta 
functions : 


where  lM  is  a function  of  {eA,  *A)  and  notation  has  been  con- 
densed by  defining 


Invoking  the  concept  of  translational  invariance  and  then 
performing  the  {R)  integrations  gives: 


l'1lr+iBM,'iA«l|  ,6‘3) 


where  r a |r-r'|.  The  resultant  quantity,  r+tBMi _tAM'  oah 
be  expressed  in  terras  of  the  known  factors  eA,  eg,  *A,  tB, 
and  r.  Elementary  geometric  considerations  lead  to: 

lr+lW  = [<r+lBM,Sin  *Bsin  eB'lAMSin  ®ASin  eA>2+ 


where  r = |r|,  *B„.  = |tBM.  |,  and  = l^l. 

Examination  of  the  key  result,  eqn.  (6-3),  indicates 
that  through  the  projection  of  molecular  to  site  correlation 
functions,  some  information  is  lost.  Thus,  given  the 
* gAMBM 1 ( r ) ) , there  is  no  unique  way  to  reconstruct  the 
9mm ' ( R ) • This  result  indicates  that  even  within  the 
concentric  shell  model,  establishing  a rigorous  means  to 
combine  site-site  correlation  functions  to  yield  molecular 
analogues  is  unlikely.  However,  numerical  calculations  may 
provide  data  from  which  approximate  relations  may  be 
developed.  In  this  effort,  simulations  are  preferable  to 
performing  the  four  dimensional  integrals  of  eqn.  (6-3). 
Computer  simulation  can  generate  the  (gra^lr)}  and  g„„,(R) 
directly.  Future  calculations  along  these  lines  is 


recommended  as  the  resulting  approximate  relations  may  lead 
to  useful  descriptions  of  dense  fluids  based  on  the  site- 


for 


this. 


were  used  to  examine  propane  data.  From 
parameter  values  were  obtained  for  CH2  groups.  This  set  of 
site  parameters  was  then  used  to  predict  PVT  data  for 
butane.  The  comparison  to  experimental  data  was  favorable, 
indicating  transferability. 

Molecular  dynamics  simulations  were  performed  to 
explore  the  feasibility  of  using  the  concentric  shell 

Lennard-Jones  simulations  which  match  liquid  nitrogen  data 
indicates  that  thermodynamic  energy  and  its  variation  with 

pressures  were  not  as  accurately  reproduced  but  this  can 
likely  be  attributed  to  larger  fluctuations  in  the  virial 
function  for  the  concentric  shell  model. 

As  an  alternative  to  simulation  for  calculating 

action  were  suggested.  These  included  perturbation  theory 
calculations  for  the  pair  correlation  function  using  the 
Lennard-Jones  fluid  as  a reference  and  using  the  concentric 
shell  model  to  develop  approximate  relations  between  site 
and  molecular  correlation  functions.  Such  relations  may 
prove  useful  in  developing  group  contribution  methods  based 
on  the  site  version  of  the  Ornstein-Zernike  equation.  Each 
of  these  alternatives  is  worthy  of  consideration  as  they  may 
lead  to  useful  correlations  for  dense  fluids. 


CHAPTER  VIZ 
CONCLUSIONS 


contribution  method  less  restrictive  than  the  currently  used 
solution  of  group  methods.  The  approach  consisted  of 
exploring  the  available  molecular  theories  that  might  serve 
as  the  basis  for  a group  contribution  method;  carefully 
analyzing  the  solution-of-groups  methods  in  order  to  cast 
the  implicit  assumptions  in  standard  statistical  mechanical 
terminology;  and  formulation  of  an  alternative  based  on  a 
concentric  shell  pair  potential  model.  The  last  consisted 
of  development  and  evaluation  of  the  potential,  derivation 
of  a group  contribution  equation  of  state  based  on  the  model 
potential,  and  molecular  dynamics  calculations  of  PVT 
properties  of  model  fluids  interacting  through  the  shell 
potential  in  both  teh  vapor  and  liquid  regions  of  the  phase 

Examination  of  rigorous  molecular  theories  such  as 
perturbation  and  integral  equation  theories  indicates  that 
their  practical  utility  is  likely  to  remain  hindered  due  to 

molecular  pair  potential  models.  Even  with  the  projected 
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The  solution -of -groups  analysis  suggests  that  some 
"smearing"  of  sites  within  molecules,  to  effectively  achieve 
site  symmetry,  can  lead  to  useful  correlations.  In  light  of 

introduced.  The  concentric  shell  model  allows  for  various 
types  of  site-site  interactions.  It  was  found  to  have 
several  attractive  properties:  the  site-site  potential  was 

recovered  in  the  limit  of  vanishing  shell  radii,-  chemically 
equivalent  site-types  can  behave  as  if  they  were  thermo- 
dynamically equivalent;  and  the  model  was  shown  to  be 

correlations  for  the  "normal"  fluids.  Thus,  one  may 

relaxed  and  the  concentric  shell  potential  forms  an  attrac- 
tive basis  for  development  of  group  contribution  methods. 

The  concentric  shell  model  was  then  utilized  in  the 
development  of  group  contribution  equations  of  state.  The 
GC-EOS  from  the  pressure  equation  was  essentially  equivalent 
to  a truncated  virial  series  unless  empirical  modifications 
are  introduced.  The  GC-EOS  of  the  van  der  Waals  form  was 
developed  more  extensively.  It  was  found  to  offer  several 
advantages  over  currently  available  GC-EOS:  an  analytic 

formula  for  the  "a"  parameter  was  developed  that  has  the 
correct  mole  fraction  dependence;  the  "b"  parameter  is  given 

need  for  mixing  rules  for  unlike  contributions  to  the 


"a"  parameter;  and  a £orm  for  the  temperature  dependence  of 
the  "a"  parameter  was  derived  from  a high  temperature 
expansion  of  the  pair  correlation  function.  These  features 
improve  upon  the  available  group  contribution  equations  of 

To  establish  credibility  for  the  model,  PVT  properties 
were  calculated  and  compared  to  experimental  data.  The 
properties  of  neopentane,  propane,  and  butane  were  generally 
well  represented  by  the  van  der  Waals  GC-EOS.  Furthermore, 
transferability  of  site  parameters  was  established.  As  is 
typical  of  equations  of  state,  some  difficulty  was 
encountered  in  accurately  reproducing  the  critical  region. 
Calculated  pressures  proved  sensitive  to  parameter  values 
and  liquid  densities.  Here,  the  parameters  used  did  not 
accurately  reproduce  the  phase  envelope.  Molecular  dynamics 
simulations  were  performed  to  explore  the  dense  fluid 

favorably  to  results  for  nitrogen,  suggesting  the  concentric 
shell  model  may  indeed  be  useful  in  modelling  dense  fluids 


Initial  evaluat 
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shell  model  were  outlined  and  are  also  worthy  of  further 
exploration. 


APPENDIX  A 

INTRAMOLECULAR  AND  PURE  COMPONENT 
CONTRIBUTIONS  TO  THE  CONFIGURATIONAL  ENERGY 


In  this  appendix  the  development  of  the  intramolecular 
expressions  for  the  internal  contributions  to  the  configura- 
tional energy  is  outlined.  Given  an  arbitrary  system 
configuration,  the  total  intramolecular  potential  energy  can 


type  M: 


ensemble  average 


°O0 ' uoM0M 1 roiM'  rBiM ' 


■ ; £ "•  J " J a 


where  the  last  equality  arises  due  to  indistinguishabillty. 
Now  define  the  intramolecular  site-site  correlation  function 
saM0M  by 

pMsaM0Mlr'  r')  S (1  - S^XH^Ir  - r^ldlr'  - rpH)> 

which  is  the  joint  probability  density  of  finding  two 
different  sites  (regardless  of  type),  o and  0,  in  the  same 
molecule  of  type  M at  r and  r' , respectively  (Chandler, 

1982).  By  definition,  the  dimensions  on  s^alr)  differ 
from  its  intermoiecular  counterpart  gaMgM . ( r ) . Then 


which  is  utilised  in  the  main  body.  For  the  pure  component 
intramolecular  potential  energy,  one  can  write,  for  type  M 


molecules : 


3aB  uaM|3M 


resulting  in  the  degenerate  form  of  (A-3),  namely, 

lirn^U  > - 2tiNm  E E | uaM3M<r)saHBM(r)r2  dr  (A-4) 

In  a similar  vein,  the  intermolecular  potential  for  the 
pure  component  case  of  type  M molecules  can  be  written 

VI  ,U  ’ " 2 i j 11  " 6i3’a  P UoMBMlroiM'rBjM* 
leading  to  the  pure  component  version  of  eqn.  (3-8): 


I uaMBM(r,9aMBM(rlr2<S 


s consistent  w 


e components,  the  derivation  here 
xture  equations  derived  separately. 


APPENDIX  B 

ALTERNATE  DERIVATION  OF 
EQUATION  (3-8) 


expression  for  the  energy,  egn.  (3-8),  can  be  obtained  by 
projecting  the  molecule-molecule  form  of  the  energy  equation 


into  site  "space".  The  essence  of  the  problem  is  to  relate 
the  molecule-molecule  pair  correlation  function,  which  is 
denoted  by,  <V  Si'1'  to  «*>  analogous  site- 


site  correlation  function,  gaMpM  (r,  r1).  This  can  be 
accomplished  by  first  defining  RM  as  the  vector  from  the 

M,  nH  as  the  set  of  (normalized)  angular  coordinates 


describing  the  external  orientation  of  the  molecule  of 
type  M,  and  l^a)  as  the  vector  from  the  center  of  mass  to 


yields 


3-2) 


■ «»»’  ('BM'  - ««  * >„<">  - V'™'  v I 

where  use  is  made  of  the  relation  for  any  two  specified 
molecules  of  type  M and  M' 


"m1  ‘ = rpM'  ‘ roM  + loM  ‘ *BM' 

Here  a is  a site  on  a molecule  of  type  M and  similarly  B is 
on  M'.  Now  average  over  external  orientations  to  obtain  the 
site-site  distribution  function: 


J »» I J I (v  % =».)  • 


% 


Here  r is  the  location  of  a,  which  is  on  M,  ana  r'  is  the 
location  of  B which  is  on  M' , The  physical  interpretation 
(B-3)  is  that  the  intersite  distance,  r - r1 , is  fixed 
via  introduction  of  the  Dirac  delta  functions,  6,  and  all 
orientations  of  M and  M1  are  averaged  INezbeda,  1977;  Gray 
and  Gubbins,  1984). 

For  an  arbitrary  configuration  of  molecules,  the  energy 


appears 


If  the  molecule-molecule  pair  potential,  u^MjM. , is  of  the 


UiMjM,[RiM'RjM''°iM'njM,J  - * | UoMpM’[roiM'rBjM'  J 


U]- 


Where  ui  and  u'  are  dummy  angles  for  orientation  averaging 
that  are  different  from  SI  and  Q'  because  the  separations 
are  r and  r'  rather  than  R and  R*.  Using  the  definition  of 
the  molecule-molecule  pair  correlation  function  (Gray  and 
Gubbins,  1984),  eqn.  (B-4)  becomes 


(B-5) 


Now  performing  the  angle  integrations,  and  noting  that  the 
distance  r'  - r is  fixed  during  these  angle  integrations, 


which  matches  egn.  (3-8).  These  derivations  from  different 
points  of  view  affirm  the  equivalence  of  the  molecular  and 
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APPENDIX  C 

SATISFYING  THE  SOLUTION- 
OF-GROUPS  EQUATIONS 

Here  it  is  shown  that  eqns.  (3-10)  and  (3-11)  will 
satisfy  eqn.  (3-1).  start  with  the  operations  prescribed  by 
equation  (3-2a)  using  eqn.  (3-9): 


•Mif' 


? J uaM 


+ lira 

V1 


-]*© 
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Substitute  egns.  (3-10)  and  (3-11)  into  eqn.  (3-9)  and  use 
egn.  (3-3)  to  obtain  AR  in  terms  of  site  types  only: 


j(c)g? 


(C-3) 


This  allows  one  to  take  the  operation  of  eqn.  ( 3— 2b) : 


[2  * * WB  f UAB(r)9AB(r)r2 


I “bc1*1 


39BCIr) 


*0 


1/Tq  B C BM  C 


[I  “BClr» 


component  limit 


jj  vamnb  1 


V1 


1 “bc1*1 


*r 

1/T0 


[r 


Usclr) 


If  one  now  views  the  distribution  functions  gM  and  s-,g  as 
functions  of  the  number  of  sites  of  various  types  (which  in 
turn  depend  on  the  number  of  molecules),  one  can  use  the 
chain  rule  in  the  form 


which  applies  to  both  and  s^g.  Equation  (2-1)  is 
satisfied  when  eqns.  (2-3)  and  (C-6)  are  used  to  show  that 
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egn.  { C— 1 ) minus  eqn.  (C-2)  is  equal  to  eqn.  ( C— 4 ) minus 
eqn.  (C-5)  when  eqns.  (2-10)  and  (2-11)  hold  and  when 


■ lim 

V1 


APPENDIX  D 

EVALUATION  OF  INTEGRALS  REQUIRED  FOR 
THE  VAN  DER  WAALS  GC-EOS 

In  this  appendix  the  integrals  of  egn.  (5-17)  are 
evaluated  using  the  shell-shell  potential  egn.  (4-19). 

Terms  of  like  powers  of  Rx  through  R4  (defined  below  egn. 
(4-18)  in  the  main  body)  will  be  considered  separately.  All 
of  the  integrals  are  of  the  same  generic  form  that  can  be 
evaluated  analytically  (Petit  Bols,  1961,-  Gradshteyn  and 
Ryzhik,  1980)  using  the  formula: 


k! (M-k) ! (p-M-l+k) 


Temporarily  neglect  all  factors  not  containing  R and 


concentrate  on  the  integrals  over  Rj  through  R4. 
Considering  in  turn  the  Rx,  r2,  R3  and  R4  terms. 


r »r« 

dMM’ 


- J " *1"“ 

dMH' 


-r  «r 

dMM' 


- r 

dMM' 


(D-5) 
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The  sum  of  egns.  (D-2)  through  ( D— 5 ) may  be  written  as 


dMMl 


where  Jmm'xi  and  Lmm'Ici  are  defined  in  eqns.  (5-19).  The 
other  contributions  to  the  potential,  and  consequently  the 
"a"  parameter,  involve  exactly  the  same  procedure  except 
that  the  value  of  p in  eqn.  ( D— 1 ) varies  with  the  type 
of  contribution  (e.g.  dipole  or  quadropole).  Terms 
involving  Rj  3 through  R4“3  integrate  to  give: 


J Rj3  R dR  = 

dim- 


f »;’»«* fi  - _1~W — .] 

(dMM'+taH+e0M'  1 [ ’J 


(D-10) 


where  eqns.  (D-8)  and  (D-9)  are  not  explicitly  shown  since 
they  differ  from  (D-6)  precisely  as  (D-3)  and  (D-4)  differ 
from  (D-2).  The  same  pattern  in  the  signs  of  the  and 
lBM  parameters  will  occur  in  all  the  integrations  as  a 
result  of  the  definitions  of  Rx  through  R4.  Again,  the  sum 
of  (D-7)  through  (D-10)  can  be  written  in  condensed  form 
with  the  aid  of  egns.  (5-19): 


f 

dMM' 


The  term  involving  R^s  integrates  as  follows 

f «;5 *«*-- — , [‘--'w-v  lit-n, 

dMM'  ,dMM1+Wt0M'>  [3  4(dMM'+laM+llPM’,J 

while  the  integrals  of  Rj5 ,R^5  and  r‘5  yield  expressions 

to  (D-12)  except  for  sign  changes  in  and  . 


identical 


The  term  of  the  shell-shell  potential  involving  r'7 
integrates  to  give: 


r 

dMM' 


while  the  Rj7,  Rj7and  R"7  terms 
fashion.  The  sum  of  these  terms 


follow  in  an  analogous 
gives  the  result 


By  including  t 
substituted  in 


proper  coefficients,  r,  from  eqn.  (4-19), 
(D-6),  (D-ll),  (D— 13 ) and  (D-15)  when 
eqn.  (5-17)  yields  (5-18)  in  the  main  body. 


(5-6)  into  eqn.  (5-26)  are 
of  eqn.  (4-19)  is  obtained: 
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Following  Che  format  of  Appendix  D,  the  remaining  terms  of 
like  order  will  be  considered  separately,  terms  not 
involving  R will  temporarily  be  suppressed,  and  eqn.  (D-l) 
will  be  used  again.  Evaluating  the  terms  in  order,  one 
finds  for  the  first  term: 


[— ‘“"‘V 

U“-S1  ln'4|ldMM'+WlBM'1 


at  which  point  it  is  convenient  to  note  that  the  symmetry  of 
the  problem  is  identical  to  that  encountered  in  Appendix  D. 
Thus  the  integrals  involving  r|-",  R*_n  and  rJ"“  will  yield 
egn.  (E-3)  except  that  the  signs  of  and  lgM.  will  change 
so  as  to  match  their  value  in  the  S.L.  The  Rx  through  R4 
integrals  can  be  summed  as  in  Appendix  D to  yield: 

0 1 ! 1 * "* 


i-iri-i)' 


The  terra  involving  R^4  integrates  as  follows: 


[l  - — V^'1  + ] 

I-  l^MM,+taM+l3M'  ^ 3(dMM,+*aM+lPM’ 1 -I 

which  can  be  summed  with  the  other  R-4  terras  to  yield: 

1 (R14‘R24"R34+R44]  *2  « = 

dMM' 


The  term  involving  Rr6  involves 
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which  can  be  summed 


r (■ 

‘W 


dR  = 


(-l)k(-l)1 


The  calculations  involving  r"8  give: 


“mM'  ,dMM,+taM+tBM')5 

r ; . '■.««■»  , 'ww2  | 
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Including  the  appropriate  coefficients  from  eqn.  (E-l)  and 
then  summing  eqns.  (E-4),  (E-5),  (E-6)  and  (E-7)  along  with 
eqn.  { E— 2 ) yield  eqn.  (5-27)  in  the  main  body. 


APPENDIX  F 

COMPUTER  PROGRAM  FOR  THE  GC-EOS 


IMPLICIT  REAL* 8 (A-H,0-Z) 

DIMENSION  V(IO) ,P(10) ,TEMP(10) ,DELH(10) 
DIMENSION  H (10) , PEXP(IO) ,HEXP(10)  DELP(IO) 
COMMON/ PARAM/S IG(2 , 2) , RNU (2)  RAD ( 2 1 RK ' EPS ( 2 2) 
COMMON/SIGMA/SIG11 , SIG22 , SIG12  ' 

OPEN ( UNIT-2 1 , FILE- ' EOS . DAT • , ACCESS- 1 SEQUENTIAL  • . 
1STATUS- 1 NEW 1 ) 

OPEN ( 22 , ACCESS- ■ SEQUENTIAL 1 , FILE- • pvt . DAT • 
1STATUS— 1 OLD 1 ) ' 

-READ  IN  TEMPERATURES,  VOLUMES,  AND  PRESSURES 
READ (22,*)  K 

REA£>(22-*)  (TEMP(I)  ,V(I)  ,PEXP(I)  ,HEXP(I)  , 1=1, K) 


C— GUESS  VALE 
123  FORMAT {' 


K-PARAMETER: 


ft  POINTS  AVAILABLE  A 


S TEMPERATURE 


GET  PARAMETERS  FOR  I 
CALL  BPARAM(D,B) 
CALL  APARAM(D.A) 
WRITE (*, 697)0 


CM** 3/MOLE) 
CALCULATE  PRESSURE  A 
T=TEMP(I) ' 


**  3/MOLECULE > 


) HELMHOLTZ  F 


HIG— RT*DL0G  ( ARG) 

CALCULATE  FREE  ENERGY  PROPER  U 
CAL/MOLE) 

HVDW=-RT*DL0G  (VOL-B)  -A/VOL 
H(I)=(HVDW-HIG)/4. 184 
IF  (HEXP(I) .EQ. 0.000)  THEN 
DELH(I)-H(I) 


DELH(I)=(H 

CONTINUE 


[)-HEXP(I))/HEXP(I) 


174 


FORMAT (8 
WRITE (21,100)  (V(I) 
FORMAT (2X, F12 . 4 ( 4X, 


P(I) , PEXP ( I ) ,H(I) ,HEXP(I) , 


WRITE(*,111) 
FORMAT (2X, ' 
WRITE ( *, 112) ( 
WRITE (21, 111) 
WRITE (2 1,112) 
FORMAT (2X, Ell 


WANNA  CONTINUE.. 


SUBROUTINE  APARAM(D.A) 


INTERACTIONS  W 


R SHELLS .... 


CT-(RNU(I)*RNU(U))/(4. 

IG=SIG(I,J) 

4R=(SSIG*SSIG)*EPS(I,J 


~ilp^ 
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APPENDIX  G 

MOLECULAR  DYNAMICS  COMPUTER  PROGRAM 


SHIFTED-FORCE  SHELL-SHELL  POTENTIAL 

SHIFTED-FORCE  CORRECTIONS  

NOT  CALCULATED. 


MOTION. 

PROPERTIES^ARE 


TRAJECTORIES  ARE  NOT 
A VERLET  NEIGHBOR- LI 


"INTVEL"  I 
:V  1/30/85)1 
GENERATOR) I 


C (THIS  VERSION 
C VARIABLES 

C TR  DESIRED  VALUE  OF  REDUCED  TEMPERATURE.  I 
C DR  REDUCED  NUMBER  DENSITY. 

C NP  NUMBER  OF  MOLECULES. 

C RC  CUTOFF  DISTANCE  FOR  POTENTIAL. 

0 DELTA  REDUCED  TIME-STEP. 

“ ixiniiiiiimiiiiiiiiiiiiiiiiiiiiiiiiiixiuiH 


INDICATES  WHETHER  EQUILIBRATION  FINISHED; 
IFLG  = 0.  EQUILIBRATION  NOT  COMPLETE. 

IFLG  = 1,  EQUILIBRATION  COMPLETED. 


1E-STEPS  PERFORMED. 
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PERFORMED.  j 

KWRITE: TIME-STEP  INCREMENT  AT  WHICH  PROPERTIES  I 


MAXKB:  TOTAL  NUMBER  OF  TIME-STEPS  TO  BE 

CALCULATED  AFTER  EQUILIBRATION.  I 

TDIST : RUNNING  VALUE  OF  MEAN-SQUARE  DISPLACEMENT. I 

XDIST : VALUE  OF  TDIST  AT  WHICH  EQUILIBRATION  IS  I 
TO  BE  TERMINATED 

ISOT : FLAG  FOR  RUNNING  EITHER  ISOTHERMALLY  (NVT)  I 
OR  AT  CONSTANT  ENERGY  (NVE) . ISOT-1 — > NVT; I 
WHILE  IF  ISOT=0  — > AN  NVE  SIMULATION.  I 


IIIIIIIIIIIIIIIIIIIIHIIIIXIXXI11IIII1IIII1I1II 
IMPLICITREAL*8  (A-H.O-Z) 

DOUBLE  PRECISION  DSQRT, DFLOAT 
COMMON/POS/XO ( 256 ) , YO (256) , ZO (256) 
COMMON/VEL/X1 (256) , Y1 (256) , Z1 (256) 
COMMON/DER/X2 (256) , Y2 (256) , Z2 (256) ,X3(256) . 

A Y3 (256) ,Z3(256),X4(256),Y4(256),Z4(256) , 

B X5(256) ,Y5(256) ,Z5(256) 

COMMON/FOR/FX(256) , FY (256) , FZ (256) 

COMMON/DISP/DAX ( 256 ) , DAY ( 25 " ' 

A X0L(256) ,Y0L(256)  - 

COMMON/NABLST/LIS 

COMMON/ PARM/FO  2 , F__ 

COMMON/PROP/IDIST(300) ,SUME,SUMV  XSUI 
COMMON/PROPA/RDEL, RMAX, RDMAX, RLIST  F' 

A ESHFT , ESHFTA, NP1 , NP2 , KSORT 

COMMON/SHELL/RAD,  SIG12 , SIG6 , RAD2 , RAD! 

COMMON/RADIAL/AVRDF ( 150 ) ,AVRRR(150) .[ 
COMMON/ISO/ISOT, IPRT 

“ SET  NUMBER  OF  PARTICLES  IN  PRIMARY  CEI 


■ SET^DESIRED  FLUID  STATE  CONDITION 

* FOR  AN  ISOTHERMAL  RUN  (NVT)  SET  ISOT=l; 


PARAMETERS 


KWRITE-10 
MAXKB-1500 
TDIST=0 . 
XDIST=0 . 4 

SET  VALUES  O 


PHYSICAL  CONSTANTS 


RADIUS3. 36 
RAD=RADIUS/SIGMA 
FACT»RNU*RNU*REPSI/  (R 


T=TR*EPSI 
VOL=PART/  DR 
VSCALE=THIRD/TR 
VOLCC=AVO* (SIGMA*I . D-8 ) **3/DR 
AND  ITS  MULTIPLES 

DELSQ-DELTA*DELTA 
DELTSQ= . 5*DELSQ 

TSTEP=DSQRT (WTMOL*SIGMA**2/AVO/EPSI/BOLZ ) 


H PREDICTOR-CORRECTOR  METHOD 


VERLET, 


R VELOCITIES  DURING  EQUIL. 


RTR=DELTA*DSQRT ( 3 . *TR) 

INCREMENT  FOR  SAMPLING  G (R)  & DIVISOR  F 
CALCULATION  OF  THE  AVERAGE  G(R) 
NY=CUBE2/RDEL-1 


REP-  (SIG12*  (R19-R29-: 
ATT=(SIG6* (R13-R23-R 
URC= ( FACT* ( REP-ATT ) ) . 
REP— (SIG12* (R110-R21< 
ATT— (SIG6* (R14-R24-R: 
DEV-FACT* (REP-ATT) 
DURC-- (URC+DEV) /RC 


-R310+R410) )/10 


CORRECTIONS  F 


S INTERACTIONS 


+(RAD2*R4)/2.) 

L 2 * ( RT1+RT2+RT3  +RT4 ) / 90 . 
>* (AT1+AT2+AT3+AT4 ) /12 . 
'2 . *PI*DR* (TERM1-TERM2 ) 


RTERM1=TERM1 
ATERM1=TERM2 
RT1=R17* (SEVEN-0. 


RTERM2=SIG12* (RT1+RT2+RT3+RT4?/10 . 
ATERM2-SIG6* (AT1+AT2+AT3+AT4 ) /4  . 
RTERM= (RTERM1+RTERM2 ) 

ATERM- (ATERM1+ATERM2 ) 


INITIALIZE  SUM  ACCUMMULATORS 


IDIST(K) *0 
PRINT  PARAMETERS 
CALL  PRNT(NP,EPSI, SIGMA, TR,T, DR, 
RMAX , RC , RDEL , RLIST , RDMAX , DELTA , N 
DE , DP , SIGMAS , EPSIS , RAD) 
IF(ISOT.EQ.l)  WRITE(6,933) 

FORMAT (T20 , 4 " •■*••• 

IF (ISOT. EQ. 1 
FORMAT (T20, ’ 

IF (IS 


4 ISOTHERMAL  R 


UN-TABLE  HEADING 
WRITE (6,930) 

FORMAT (1H1///7X,  'KB',2X,  'TIME', 
■PI' ,5X, 'ENRG ' , 6X, 'El' ,7X, 'DISI 
4X, 'NAB' ,3X, 'ETOTAL'/) 

^3AD  INITIAL  POSITIONS  OF  ATOMS 


LOAD  INITIAL  VELOCITIES  OF  ATO 
CALL  INTVEL ( AHEAT , RTR , PART , N 
ASSIGN  INITIAL  ACCELERATIONS  B 


SCALE  ACCELERATIONS  & 

X2(I)=FX(I)*DELTSQ 
Y2(I)=FY(I)*DELTSQ 
Z2 ( I ) =FZ ( I ) * DELTSQ 


E STARTING  POSITIONS 


CONTINUE 

NSTART-1 
DO  599  NTIMES=NSTART,MAXKB 


F SIMULATION 


N SQUARE  DISPLACEMENT  & 


TDIST=TDIST+DAX(I) **2+DAY(I) **2+DAZ ( 

SUMVEL=SUMVEL+X1(I) **2+Yl{I) **2+Zl!I 

CONTINUE 

TDIST-TDIST/PART 

"K=SUMVEL/ ( 2 . *PART*DELSQ) 

3 PROPERTY  AVERAGES 

M+SUMVEL 


ACCUMMULATE  S 


SUMV=SUMV+TOTV 

PROPERTY  CALCULATION  & PRINT-OUT  AT  INTERVAL! 
IF(MOD(KB,KWRITE) .NE.O)  GOTO  550 
FKB=DFLOAT(KB) *PART 
TMP=XSUM/ (3 . *DELSQ*FKB) 

ENR=  ( SUME/  FKB+CORE ) 

VIR=VSCALE*  ( SUMV/FKB+CORV) 

PRES=TMP*DR* ( 1 . -VIR) 

El-TOTE/PART+CORE 

ET0T=E1+EK 

P1=TMP*DR*  ( 1 . -VSCALE*  (TOTV/PART+CORV)  ) 
RLTIM=DELTA*DFLOAT(KB) *TSTEP 
WRITE (6,940)  KB, RLTIM , PRES , PI , ENR , El , 
TDIST , TMP , NABTOT , ETOT 

FORMAT ( 1H  ,3X,I5.F7.3.5F9  3 F9  4 TB  IV 


4.9) 


G(R)  i 


' INTERVALS 


IF (MOD (KB, 500) . EQ. O.OR. KB. EQ.MAXKB) 

A CALL  RDF(DR,TMP,RDEL,NP,NY,KB) 

CONTINUE 

DURING ^FIRST^OF  RUN,  SCALE  VELOCITIES  FOR  TE 
A IFLG  n^kb)'  (SUMVEL,  AH  EAT , TDIST , XDIST , NP , 
L EQBRAT  ( SUMVEL,  AHEAT , 


TDIST, XC 
CONTINUE 


?,IFLG,N 


) RRR(I) , AVRDF ( I ) 


NSTART-KB+I 

SUBROUTINE  FCC(CUBE,NP) 
CALCULATE  POSITIONS  OF  SITES  01 
CUBIC  LATTICE 

BASED  ON  A BOX  OF  SIDE  = 100 
IMPLICITREAL*8  (A-H,0-Z) 
COMMON/ POS/XO (256) ,Y0(256) ,2( 
NC=(DFLOAT(NP)/4.)**(l./3. )+. 
XL=100./DFLOAT(NC) 


<0(l)-0. 

fO(l)=0. 


SUBROUTINE  INTVEL  ( AHEAT , R 
ASSIGN  INITIAL  VELOCITIES  T 
IMPLICITREAL*8  (A-H,0-Z) 
DOUBLE  PRECISION  DSQRT.DF 
COMMON/VEL/X1 (256) ,Y1(256 
MS EE 0=34567 


S THE  NERDC  RANDOM  # GENERATOR; 


0=RNDMF(B) 
A+. 74864321 
LL  RSEED(A) 


2 RANDOM  i GENERATOR 


T TOTAL  MOMENTUM=ZE 


XI ( I ) =X1 ( I ) -SUMX/PART 
Y1 (I) -Y1 ( I) -SUMY/PART 
Z1(I)=Z1(I) -SUMZ/PART 


SUBROUTINE  PREDCT(NP) 

USE  FIFTH-ORDER  TAYLOR  SERIES  T 
POSITIONS  S THEIR  DERIVATIVES  A 
IMPLICITREAL*8  (A-H.O-Z) 
COMMON/POS/XO (2S6) , YO (256) ,Z0( 
COMMON/VEL/X1 (256) ,Y1(256) ,Z1( 
COMMON/DER/X2 (256) ,Y2(256) ,Z2( 
Y3 (256) , Z3 (256) ,X4(256) ,Y4(256 
# Y5 (256) , Z5 (256) 

COMMON/ FOR/ FX ( 2 5 6 ) ,FY{256) ,FZ( 


C3(I)+X4(I)+X5(I) 


!0(I)+Z1(I^+Z2(I)+Z3(D+2 


tl(I)-Yl(I)+2. 

;.*Y5(i) 

!l(I)=Zl(I)+2. 

(2(I)=X2(I)+3. 

(2(I)-Y2(I)+3. 

!2(I)-Z2(I)+3. 

C3(I)=X3(I)+4. 

f3(I>=Y3(I)+4. 


INITIALIZE  FORCE  ARRAYS 
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SUBROUTINE  EVAL(TOTV,  TOTE,  CUBE,  CUBE2  KB  NP 
A NABTOT)  ' ’ ' ' 

EVALUATE  FORCES  ON  ATOMS  WITH  PAIRWISE  ADDITIVE 

DOUBLE  PRECISION  DABS, DSQRT, DSIGN 
COMMON/POS/XO (256) , YO (256) ,Z0{256) 

COMMON/ FOR/ FX ( 2 56 ) , FY (256) , FZ (256) 
COMMON/NABLST/LIST (10000) , NABORS (256) 

COMMON/ PROP/IDIST( 300) , SUME, SUMV, XSUM  ISUM 
COMMON/PROPA/RDEL, RMAX, RDMAX, RLIST, FSHFT 
A ESHFT , ESHFTA , NP1 , NP2 , KSORT 

COMMON/SHELL/RAD, SIG12 , SIG6 , RAD2 , RAD22 , RAD23 , 


IF(MOD(KB, KSORT) .E 


E NEIGHBOR-LIST 

JBEGIN-NABORS ( I ) 
JEND  =NABORS(I+l)-; 
IF (JBEGIN. GT. JEND)  C 


STORE  POSITION  0 


IF (XR. GT. CUBE2 
IF(YR.GT.CUBE2 
IF(ZR.GT.CUBE2)  Z=i 


(XR-CUBE) * DSIGN ( 1 . D 
'‘"1-CUBE)  ‘DSIGN  ( 1 . D 
t-CUBE)  ‘DSIGN (l.D 


CALCULATE  R(IJ) 

RCTR=DSQRT(RSQ) 
INCREMENT  COUNTER  FOR  G(R' 
AT  INTERVALS 

IF(LCHK.NE.l)  GOTO  i 
IF (RSQ . GT . RDMAX ) C 
IJ-RCTR/RDEL+.5 
IDIST(IJ)=IDIST(Iu 
IF(RSQ.GT.RLIST)  C 

LIST(K)=J 

IF(RSQ.GT.RMAX)  GOTO 
EVALUATE  FORCES  ON  ATOMS 
RSI-1. /RSQ 


B^rn^R110_R210_R:310+R4l0)/lo! 
RTERM—SIG12* (RT1+RT2)  J/iu' 

£S:ffi2:5?2:2??:5i?wi2-*R> 


ACCUMMULATE  ENERGY  AND  V 
UR=(SIG12*RT1-SIG6»AT1 
TOTE-TOTE+UR+ESHFT+RCT 
TOTV=TOTV-RPL 
CONTINUE 
CONTINUE 

IF(LCHK.NE.l)  RETURN 
NABORS  (NP)=K-H 


SUBROUTINE  CORR ( DELTSQ , CUBE , NP ) 

CORRECT  PREDICTED  POSITIONS  AND  DERIVATI 
IMPLICITREAL*8  (A-H.O-Z) 

COMMON/POS/XO (256) ,Y0(256) ,Z0(256) 
COMMON/VEL/X1 (256) , Y1 (256) ,Z1(256) 
COMMON/DER/X2 (256) ,Y2(256) ,Z2(256) ,X3(2 
A Y3(256) ,Z3(256) , X4 (256) , Y4 (256) 

B ,X5(256) ,Y5(256) ,Z5(256) 

COMMON/ FOR/ FX (256) ,FY(256) ,FZ(256) 
COMMON/DISP/DAX (256) , DAY (256) , DAZ (256) 

A X0L(256) , Y0L(256) ,Z0L(256)  ' 

COMMON/PARM/F02 , F12 , F32 , F42 , F52 
XERR-X2 (I) -DELTSQ*FX(I) 

YERR=Y2 (I) -DELTSQ*FY (I) 

ZERR=Z2 (I) -DELTSQ*FZ ( I) 


DISPLACEMENTS 

DAX(I)-DAX(I)-XO 
DAY (I) =DAY (I) -Y0 
DAZ (I) -DAZ (I) -Z0 


CONSERVATION 


SUBROUTINE  EQBRAT (SUMVEL, AHEAT, TDIST, XC 


COMMON/PROP/IDIST (3C 
COHMON/ISO/ISOT,  IPRT 
IF(ISOT. EQ. 1) GOTO  721 
IF (TDIST. GT.XDIST)  GOTO  750 
IF ( TDIST . GT . XDIST . AND . IPRT . EQ . 0 ) 
HEAT=DSQRT  ( AHEAT/SUMVEL) 


F EQUILIBRATION  S 


, SET  PROPERTY 


r.100)  RETURN 


WRITE (6, 777) 
FORMAT (/10X, 
IPRT=1 

WRITE (6,930) 
FORMAT ( 1H1///7Xj 


EQUILIBRATION  COMPLETE  '//) 


SUBROUTINE  RDF ( DR , TMP , RDEL, NP  NY  KB) 
NORMALIZE  COUNTERS  FOR  RADIAL  DISt!  FUNCTI 
IMPLICITREAL*8  (A-H.O-Z) 
COMMON/PROP/IDIST(300) ,SUME,SUMV  XSUM  IS 
COMMON/RADIAL/AVRDF ( 150 ) ,AVRRR(150)  RDFN 
PRINT  HEADING 

FORMAT(1H1////T20,  'RADIAL  DISTRIBUTION 


IF(IDIST(J) .EQ.O)  GOTO  780 

V= (4 . *3 . 14159/3 . ) * ( 3 . *DFLOAT (; 
X=IDIST(J) 

GR-X/V/Y/PDEN 

RRR«RDEL»DFLOAT(J) 

"t*"‘  J,RRR,IDIST(J)  ,GF 


FORMAT(28X,I3 
IP(I . GT. 149) 
IF(RRR.LT.AVR 


CONTINUE 
WRITE(6,970) 
FORMAT (1H1///) 


192 


50)  , AVRRR  (150)  , RE 
i,X5,V5, 25/2304*0 
i , FZ/1536*0 . DO/ 


WRITE (6, 302) 
FORMAT(T20, • 
EPSILON/K  = ' 


FORMAT  (T20,  ' 
WRITE (6, 904) 
WRITE(6,922) 
FORMAT (T20, ' 


WRITE (6  ^ 926) 
FORMAT (T20, 1 


D SHELL  PARAMETERS: 


, 'ENERGY  CORRECTION  = 
, 'PRESSURE  CORRECTION  - 


F OUTPUT  COLUMN  HEADINGS 


WRITE (6,904) 

WRITE(6,902) 

PRINT  IDENTIFICATION  0 
WRITE (6, 951) 

FORMAT (/////T20,  'OUTPUT  COLUMN  HEADING 
SYMBOLS'/) 

WRITE (6,952) 

FORMAT  (T2  4 , "I 


WRITE (6, 954) 
FORMAT (T2 4, 'P! 

FORMAT (T24,  'El 
ENERGY ' ) 
WRITE(6,956) 
FORMAT (T24, 'E: 


FORMAT  (T2  4 , "I 
FORMAT (T24, 'N 


ME-ELAPSED  TIME  IN  PICOSECONDS') 

1ES— RUNNING  AVE.  PRESSURE,  P* ' ) 
-INSTANTANEOUS  PRESSURE ' ) 

RG-RUNNING  AVE  CONF  INTERNAL 

-INSTANTANEOUS  CONF  INTERNAL  ENERGY  ' ) 
ST-MEAN  SQ  DISPLACEMENT  FROM 

MP-RUNNING  AVERAGE  TEMPERATURE ' ) 
B-INSTANTANEOUS  TOTAL  NEIGHBORS 


960  RETURN <T24' 'ETOTAL-INSTANTANEOOS  TOTAL  ENERGY') 
FUNCTION  RANDOM (MSEED) 

£ SL""™  CREATES  A PSEUDO-RANDOM  NUMBER  BY  A 
n 2S»^LiCATIVE  CONGRUENTIAL  METHOD.  THE  MAXIMUM 
C PERIOD  OF  THIS  GENERATOR  IS  THE  VALUE  OF  THE  MODULUS 
C BE  1658  ™A«  » A“  DEPENDS 

c UPON  THE  MULTIPLIER  MULT.  SINCE  HERE  WE  NEED  AT  MOST 
C A FEW  HUNDRED  RANDOM  NUMBERS,  WE  HAVE  SET  M TO  i wit, 
C VALUE  AND  HAVE  CHOSEN  AN  ARBITRARY  VALUE  FOR  MULT  THE 

c ?™LcomtoterEhamSa^KES  seI  °l 

C COMPUTER^ PROGRAMMING , " VOL.  2,  ADDISON-WESLEY,  READING, 
IMPLICIT  REAL*8 * (A-H.O-2) 


MNEW  =MOD(MPROD,M) 

’AIN  A REAL  RANDOM  NUMBER  ON  THI 
■ DFLOAT ( MNEW ) / DFLOAT ( M ) 
THE  INTERVAL  -• 


RANDOM- 

MSEED-MNEW 
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